A few days ago I coded up a demo of anomaly detection using principal component analysis (PCA) reconstruction error. I implemented the PCA functionality — computation of the transformed data, the principal components, and the variance explained by each component — from semi-scratch, meaning I used the NumPy linalg (linear algebra) library eig() function to compute eigenvalues and eigenvectors.

And it was good.

But in the back of my mind, I was thinking that I should have verified my semi-from-scratch implementation of PCA because PCA is very, very complex and I could have made a mistake.

The from-scratch version (left) and the scikit version (right) are identical except that some of the transformed vectors and principal components differ by a factor of -1. This doesn’t affect anything.

So I took my original from-scratch PCA anomaly detection program and swapped out the PCA implementation from the scikit sklearn.decomposition library. And as expected, the results of the scikit-based PCA program were identical to the results of the from-scratch PCA program. Almost.

My from-scratch code looks like:

import numpy as np

def my_pca(X):
  # returns transformed X, prin components, var explained
  dim = len(X[0])  # n_cols
  means = np.mean(X, axis=0)
  z = X - means  # avoid changing X
  square_m = np.dot(z.T, z)
  (evals, evecs) = np.linalg.eig(square_m) 
  trans_x = np.dot(z, evecs[:,0:dim]) 
  prin_comp = evecs.T  
  v = np.var(trans_x, axis=0, ddof=1) 
  sv = np.sum(v)
  ve = v / sv
  # order everything based on variance explained
  ordering = np.argsort(ve)[::-1]  # sort order high to low
  trans_x = trans_x[:,ordering]
  prin_comp = prin_comp[ordering,:]
  ve = ve[ordering]
  return (trans_x, prin_comp, ve)

X = (load data from somewhere)
(trans_x, p_comp, ve) = my_pca(X)

The scikit-based code looks like:

import numpy as np
import sklearn.decomposition

X = (load data from somewhere)
pca = sklearn.decomposition.PCA().fit(X)
trans_x = pca.transform(X)
p_comp = pca.components_
ve = pca.explained_variance_ratio_

All the results were identical except that the internal transformed X values and the principal components, sometimes differed by a factor of -1. As it turns out this is OK because PCA computes variances and the sign doesn’t affect variance.

The advantage of using scikit PCA is simplicity. The advantages of using PCA from scratch are 1.) you get fine-tuned control, 2.) you remove an external dependency, 3.) you aren’t using a mysterious black box.

PCA is interesting and sometimes useful, but for tasks like dimensionality reduction and reconstruction, deep neural techniques have largely replaced PCA.

PCA was developed in 1901 by famous statistician Karl Pearson. I wonder if statisticians of that era imagined today’s deep neural technologies. Three images from the movie “Things to Come” (1936) based on the novel of the same name by author H. G. Wells.

Demo code:

# pca_recon_skikit.py
# exactly replicates iris_pca_recon.py scratch version

import numpy as np
import sklearn.decomposition

def reconstructed(X, n_comp, trans_x, p_comp):
  means = np.mean(X, axis=0)
  result = np.dot(trans_x[:,0:n_comp], p_comp[0:n_comp,:])
  result += means
  return result

def recon_error(X, XX):
  diff = X - XX
  diff_sq = diff * diff
  errs = np.sum(diff_sq, axis=1)
  return errs

def main():
  print("\nBegin Iris PCA reconstruction using scikit ")
  np.set_printoptions(formatter={'float': '{: 0.1f}'.format})

  X = np.array([
    [5.1, 3.5, 1.4, 0.2],
    [5.4, 3.9, 1.7, 0.4],

    [6.4, 3.2, 4.5, 1.5],
    [5.7, 2.8, 4.5, 1.3],

    [7.2, 3.6, 6.1, 2.5],
    [6.9, 3.2, 5.7, 2.3]])

  print("\nSource X: ")
  print(X)

  print("\nPerforming PCA computations ")
  pca = sklearn.decomposition.PCA().fit(X)
  trans_x = pca.transform(X)
  p_comp = pca.components_
  ve = pca.explained_variance_ratio_
  print("Done ")

  print("\nTransformed X: ")
  np.set_printoptions(formatter={'float': '{: 0.4f}'.format})
  print(trans_x)

  print("\nPrincipal components: ")
  np.set_printoptions(formatter={'float': '{: 0.4f}'.format})
  print(p_comp)

  print("\nVariance explained: ")
  np.set_printoptions(formatter={'float': '{: 0.5f}'.format}) 
  print(ve)

  XX = reconstructed(X, 4, trans_x, p_comp)
  print("\nReconstructed X using all components: ")
  np.set_printoptions(formatter={'float': '{: 0.2f}'.format})
  print(XX)

  XX = reconstructed(X, 1, trans_x, p_comp)
  print("\nReconstructed X using one component: ")
  np.set_printoptions(formatter={'float': '{: 0.2f}'.format})
  print(XX)

  re = recon_error(X, XX)
  print("\nReconstruction errors using one component: ")
  np.set_printoptions(formatter={'float': '{: 0.3f}'.format})
  print(re)

  print("\nEnd PCA scikit ")

if __name__ == "__main__":
  main()