ENGO 629: Robust PCA vs Standard PCA with presence of outliers

Example_Robust_PCA_with_MAD.py

#!/usr/bin/env python3

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

outliers  = np.loadtxt("outliers.txt")
plane_pts = np.loadtxt("plane_pts.txt")

plane_with_outliers = np.vstack((plane_pts, outliers))

# plot figure
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(plane_pts[:,0], plane_pts[:,1], plane_pts[:,2], 'b.')
ax.plot(outliers[:,0], outliers[:,1], outliers[:,2], 'r.')
ax.xaxis.set_label("X")
ax.yaxis.set_label("Y")
ax.zaxis.set_label("Z")
plt.show()

def cov_MAD(data):
    mu = np.median(data, axis=0)
    CV = np.zeros((data.shape[1], data.shape[1]))

    for i in range(data.shape[1]):
        for j in range(data.shape[1]):
            CV[i,j] = 1.4826 * np.median(np.abs(data[:,i] - mu[j]))
    return CV

# Print CV-matrix and plane normal using LS-estimator on non-outliers
# Note that the plane normal corresponds to the smalles eigenvalue / vector pair
# Also note that numpy is weird in that it calculates covariance row-wise,
# whereas most literature does this column-wise. Hence the .T suffix to transpose
# the matrix
print(np.cov(plane_pts.T))
print(np.linalg.eigh(np.cov(plane_pts.T)))

# Do the same for the contaminated points
print(np.cov(plane_with_outliers.T))
print(np.linalg.eigh(np.cov(plane_with_outliers.T)))

# Do the same using cov_MAD instead of LS covariance estimate
print(cov_MAD(plane_with_outliers))
print(np.linalg.eigh(cov_MAD(plane_with_outliers)))

outliers.txt

3.882364116690339984e+01 1.273181953482825968e+01 5.456523248642312041e+01
3.660421700355402663e+01 6.730261239412651264e+01 1.090938175428481998e+02
3.157262643892146770e+01 6.247625175616308013e+01 6.575782577853652811e+01
1.069563129517350220e+01 8.551648133394620288e+01 1.101410896100223624e+02
8.125926007570463128e+01 8.222708675710940440e+01 1.034688985939998815e+02
8.879870551595662675e+01 5.796213895908790903e+01 6.086179126047318277e+01
8.445116656527771681e+01 2.700858566177652875e+01 7.305038025971492743e+01
6.036967689329919295e+01 5.417916402549300159e+01 6.089150409165332434e+01
3.774255768986154180e+01 1.601301489915656262e+01 1.265015380330221291e+02
9.833126795169448542e+01 3.369413889040184795e+00 8.763290235181293042e+01
4.389778156386009300e+01 6.814961811722487539e+01 1.080043267163039076e+02
2.025713240673908544e+01 3.306349392430040268e+01 1.431349647210099647e+02
7.180638870352683512e+01 4.831031429593588200e+01 1.063373879393368213e+02

plane_pts.txt

4.783026353607691306e+01 4.222892259499848677e+01 0.000000000000000000e+00
2.747536655891493851e+01 6.274445921048194919e+01 0.000000000000000000e+00
8.960142389868082091e+01 1.509554474823211478e+01 0.000000000000000000e+00
7.056119786260236992e+01 3.332293959457660293e+01 0.000000000000000000e+00
5.874618635418080714e+01 2.935473765713727445e+01 0.000000000000000000e+00
5.331315415769025634e+01 9.315365122353155414e+01 0.000000000000000000e+00
4.412395284134726126e+01 9.591393955923560100e+01 0.000000000000000000e+00
8.835108666640995523e+01 4.699679985131254512e+01 0.000000000000000000e+00
1.352551139452153883e+01 6.281086005001023764e+00 0.000000000000000000e+00
8.243992708612790921e+01 8.599430521844130482e+00 0.000000000000000000e+00
1.012421054838408274e+00 1.298747115510776595e+01 0.000000000000000000e+00
6.268717024403708393e+01 8.500210626388889068e+01 0.000000000000000000e+00
7.351008930538954189e+01 2.766346954341744180e+01 0.000000000000000000e+00
8.275808827650590160e+01 9.840985035177662610e+01 0.000000000000000000e+00
4.867363431568371368e+01 9.403353489752504402e+01 0.000000000000000000e+00
9.035277518070317093e-01 5.815646856815433097e+01 0.000000000000000000e+00
2.487065359239920426e+01 5.022888739055871099e+01 0.000000000000000000e+00
1.082446397745904321e+01 7.863371739732087917e+01 0.000000000000000000e+00
9.267694704750030610e+01 1.364632484860770489e+01 0.000000000000000000e+00
9.733531542411522075e-01 3.234814166044264283e+01 0.000000000000000000e+00
2.220347695209803263e+01 3.593832123328961359e+01 0.000000000000000000e+00
8.681653152227688963e+01 7.039148865697329427e+01 0.000000000000000000e+00
9.131825215321454436e+01 3.655463849711681235e+01 0.000000000000000000e+00
7.716183638465525974e+01 1.840818407659181588e+01 0.000000000000000000e+00
7.978335462448590931e+01 9.972683235729597584e+01 0.000000000000000000e+00
2.732083810661735157e+01 1.468049183067124908e+01 0.000000000000000000e+00
6.763985884736212206e+00 5.003775384438180396e+01 0.000000000000000000e+00
8.370032125353708352e+01 8.226453731348799181e+01 0.000000000000000000e+00
1.067739202447044455e+01 4.566798164489856049e+01 0.000000000000000000e+00
8.813934856999232181e+01 6.311344119262513885e+01 0.000000000000000000e+00
6.098612091843565963e+01 9.956004075335748382e+01 0.000000000000000000e+00
7.044855800029642978e+01 3.926494814218253282e+01 0.000000000000000000e+00
9.345759631474568607e+01 4.772969908486979662e+00 0.000000000000000000e+00
1.895640795082850261e+01 8.407256554433337215e+01 0.000000000000000000e+00
5.936206046435822259e+01 6.956804317857641706e+01 0.000000000000000000e+00
4.400441514604012383e+01 6.896471661992332258e+01 0.000000000000000000e+00
2.219362082370686196e+01 8.462120883715090258e+01 0.000000000000000000e+00
2.603314877631357405e+01 3.295073461574082785e+01 0.000000000000000000e+00
5.907861837159400409e+01 3.813330695777558077e+01 0.000000000000000000e+00
3.761723860773469319e+00 8.405606672394218037e+01 0.000000000000000000e+00
5.061728623739494282e+01 5.883080621701250834e+01 0.000000000000000000e+00
5.237864965485015034e+01 4.344274848196572059e+00 0.000000000000000000e+00
8.554848866610279856e+01 6.798231424947989865e+01 0.000000000000000000e+00
9.126098442603142757e+01 9.095736275155907435e+01 0.000000000000000000e+00
7.159528792399333952e+01 6.013092435454478135e+01 0.000000000000000000e+00
1.214927166604680053e+01 8.325346570351666031e+00 0.000000000000000000e+00
9.985656909807927306e+00 9.977030667586731738e+01 0.000000000000000000e+00
7.205582491863204986e-01 1.423206104407871386e+01 0.000000000000000000e+00
6.684416700741560646e+01 6.602683431069176834e+01 0.000000000000000000e+00
7.366385436488089056e+01 9.220175623852912850e+01 0.000000000000000000e+00
8.154599109535753954e-01 2.152131521233397038e+01 0.000000000000000000e+00
2.449986961416599129e+01 8.346758678736911463e+01 0.000000000000000000e+00
8.977285531576903566e+01 1.889206000936143326e+01 0.000000000000000000e+00
4.378721179030906541e+01 5.870662565512819953e+01 0.000000000000000000e+00
2.253083652873287690e+01 9.006200742093161793e+00 0.000000000000000000e+00
3.245064825601306779e+01 3.608392389176800918e+01 0.000000000000000000e+00
8.476120845672376447e+01 2.630335494200716084e+01 0.000000000000000000e+00
7.197195679250922851e-01 1.264577383054389337e+01 0.000000000000000000e+00
3.985885526764296571e+01 7.745070287923883257e+01 0.000000000000000000e+00
5.615296374588199058e+01 8.415449440547763515e+01 0.000000000000000000e+00
3.488296704646150914e+01 6.035161388360506152e+01 0.000000000000000000e+00
4.682125390091221107e+01 8.588927439201225411e+01 0.000000000000000000e+00
7.688848232609009870e+01 1.875982285124633275e+00 0.000000000000000000e+00
9.217829444134140715e+01 1.885703159475632518e+01 0.000000000000000000e+00
3.443779849412486271e+01 5.989541437417845771e+01 0.000000000000000000e+00
1.644977797757662952e+00 7.946575475804479538e+01 0.000000000000000000e+00
1.934446643721993198e+01 8.044842491280380159e+01 0.000000000000000000e+00
1.681025567627382244e+01 1.688947430292439122e+01 0.000000000000000000e+00
9.611793622003700932e+01 7.213400548610461271e+01 0.000000000000000000e+00
8.213551201910430422e+00 7.112020666946095560e+01 0.000000000000000000e+00
7.378107502225100234e+01 6.170587784477304893e+01 0.000000000000000000e+00
7.363833862359376781e+01 3.709210599380254081e+01 0.000000000000000000e+00
8.575022209175347143e+01 7.117222324602508365e+01 0.000000000000000000e+00
2.346789057994688221e+01 1.640377174568887853e+01 0.000000000000000000e+00
1.441905999800263771e+01 6.137083395227771376e+01 0.000000000000000000e+00
4.823536853621803999e+01 6.367816817765719151e+01 0.000000000000000000e+00
7.759488548277620623e+01 8.223660540984456091e+00 0.000000000000000000e+00
1.887924527365631988e+01 9.677486792108517477e+01 0.000000000000000000e+00
8.930171998973163738e+01 9.964680276150726002e+01 0.000000000000000000e+00
4.819926488756603078e+01 7.271098486082938450e+01 0.000000000000000000e+00
5.908768321571532312e+01 9.937923702698313377e+01 0.000000000000000000e+00
8.940660626660165633e+01 6.761115992621826365e+01 0.000000000000000000e+00
5.526108809939371014e+01 8.095858968258502841e+01 0.000000000000000000e+00
4.654840828267331432e+01 7.049801724465754660e+01 0.000000000000000000e+00
2.051424265704519101e+01 1.089074668189141626e+01 0.000000000000000000e+00
5.618756155021387144e-01 9.822285118193298104e+01 0.000000000000000000e+00
9.373290067109671497e+01 5.725632358352798690e+01 0.000000000000000000e+00
1.619451958403292835e+01 4.350516438830279498e+01 0.000000000000000000e+00
9.935960357577704372e+01 9.588615527820198281e+01 0.000000000000000000e+00
3.215981637362123990e+01 3.782781211761505347e+01 0.000000000000000000e+00
4.314143823532593558e+01 6.820926032245073145e+01 0.000000000000000000e+00
6.110631090713187774e+01 1.640461839854285842e+01 0.000000000000000000e+00
2.527330437377221983e+01 6.498053417383609087e+01 0.000000000000000000e+00
6.884643875145147263e+01 8.691486455201844308e+01 0.000000000000000000e+00
1.622052030225761143e+01 2.938623325201211856e+01 0.000000000000000000e+00
5.842776745658516546e+01 3.340096636669720453e+00 0.000000000000000000e+00
2.706265692915742349e+01 1.833781674865477740e+01 0.000000000000000000e+00
6.395455379675954077e+01 7.912679589460034890e+01 0.000000000000000000e+00
1.593782525563480768e+01 2.806652393128061718e+01 0.000000000000000000e+00
9.199081246328051975e+01 8.040742650052095541e+01 0.000000000000000000e+00