Example Principal Component Analysis - Revision history

Peter at 15:32, 7 July 2021

2021-07-07T15:32:45Z

Show changes

Peter at 19:45, 18 May 2017

2017-05-18T19:45:21Z

Peter at 00:59, 10 August 2015

2015-08-10T00:59:18Z

Peter at 02:41, 1 April 2015

2015-04-01T02:41:23Z

Peter: Created page with "Principal Component Analysis (PCA) is a popular and simple to implement classification technique, often used in face recognition. The following is an example of how to implem..."

2015-03-22T05:30:30Z

Created page with "Principal Component Analysis (PCA) is a popular and simple to implement classification technique, often used in face recognition. The following is an example of how to implem..."

New page

Principal Component Analysis (PCA) is a popular and simple to implement classification technique, often used in face recognition. The following is an example of how to implement it in EJML using the procedural interface. It is assumed that the reader is already familiar with PCA.

Example on GitHub:
* [https://github.com/lessthanoptimal/ejml/blob/master/examples/src/org/ejml/example/PrincipalComponentAnalysis.java PrincipalComponentAnalysis]

For additional information on PCA:
* [http://en.wikipedia.org/wiki/Principal_component_analysis General information on Wikipedia]

= Sample Code =

<syntaxhighlight lang="java">
/**
* 
* The following is a simple example of how to perform basic principal component analysis in EJML.
* 
*
* 
* Principal Component Analysis (PCA) is typically used to develop a linear model for a set of data
* (e.g. face images) which can then be used to test for membership. PCA works by converting the
* set of data to a new basis that is a subspace of the original set. The subspace is selected
* to maximize information.
* 
* 
* PCA is typically derived as an eigenvalue problem. However in this implementation {@link org.ejml.interfaces.decomposition.SingularValueDecomposition SVD}
* is used instead because it will produce a more numerically stable solution. Computation using EVD requires explicitly
* computing the variance of each sample set. The variance is computed by squaring the residual, which can
* cause loss of precision.
* 
*
* 
* Usage: 
* 1) call setup() 
* 2) For each sample (e.g. an image ) call addSample() 
* 3) After all the samples have been added call computeBasis() 
* 4) Call sampleToEigenSpace() , eigenToSampleSpace() , errorMembership() , response()
* 
*
* @author Peter Abeles
*/
public class PrincipalComponentAnalysis {

// principal component subspace is stored in the rows
private DenseMatrix64F V_t;

// how many principal components are used
private int numComponents;

// where the data is stored
private DenseMatrix64F A = new DenseMatrix64F(1,1);
private int sampleIndex;

// mean values of each element across all the samples
double mean[];

public PrincipalComponentAnalysis() {
}

/**
* Must be called before any other functions. Declares and sets up internal data structures.
*
* @param numSamples Number of samples that will be processed.
* @param sampleSize Number of elements in each sample.
*/
public void setup( int numSamples , int sampleSize ) {
mean = new double[ sampleSize ];
A.reshape(numSamples,sampleSize,false);
sampleIndex = 0;
numComponents = -1;
}

/**
* Adds a new sample of the raw data to internal data structure for later processing. All the samples
* must be added before computeBasis is called.
*
* @param sampleData Sample from original raw data.
*/
public void addSample( double[] sampleData ) {
if( A.getNumCols() != sampleData.length )
throw new IllegalArgumentException("Unexpected sample size");
if( sampleIndex >= A.getNumRows() )
throw new IllegalArgumentException("Too many samples");

for( int i = 0; i < sampleData.length; i++ ) {
A.set(sampleIndex,i,sampleData[i]);
}
sampleIndex++;
}

/**
* Computes a basis (the principal components) from the most dominant eigenvectors.
*
* @param numComponents Number of vectors it will use to describe the data. Typically much
* smaller than the number of elements in the input vector.
*/
public void computeBasis( int numComponents ) {
if( numComponents > A.getNumCols() )
throw new IllegalArgumentException("More components requested that the data's length.");
if( sampleIndex != A.getNumRows() )
throw new IllegalArgumentException("Not all the data has been added");
if( numComponents > sampleIndex )
throw new IllegalArgumentException("More data needed to compute the desired number of components");

this.numComponents = numComponents;

// compute the mean of all the samples
for( int i = 0; i < A.getNumRows(); i++ ) {
for( int j = 0; j < mean.length; j++ ) {
mean[j] += A.get(i,j);
}
}
for( int j = 0; j < mean.length; j++ ) {
mean[j] /= A.getNumRows();
}

// subtract the mean from the original data
for( int i = 0; i < A.getNumRows(); i++ ) {
for( int j = 0; j < mean.length; j++ ) {
A.set(i,j,A.get(i,j)-mean[j]);
}
}

// Compute SVD and save time by not computing U
SingularValueDecomposition<DenseMatrix64F> svd =
DecompositionFactory.svd(A.numRows, A.numCols, false, true, false);
if( !svd.decompose(A) )
throw new RuntimeException("SVD failed");

V_t = svd.getV(null,true);
DenseMatrix64F W = svd.getW(null);

// Singular values are in an arbitrary order initially
SingularOps.descendingOrder(null,false,W,V_t,true);

// strip off unneeded components and find the basis
V_t.reshape(numComponents,mean.length,true);
}

/**
* Returns a vector from the PCA's basis.
*
* @param which Which component's vector is to be returned.
* @return Vector from the PCA basis.
*/
public double[] getBasisVector( int which ) {
if( which < 0 || which >= numComponents )
throw new IllegalArgumentException("Invalid component");

DenseMatrix64F v = new DenseMatrix64F(1,A.numCols);
CommonOps.extract(V_t,which,which+1,0,A.numCols,v,0,0);

return v.data;
}

/**
* Converts a vector from sample space into eigen space.
*
* @param sampleData Sample space data.
* @return Eigen space projection.
*/
public double[] sampleToEigenSpace( double[] sampleData ) {
if( sampleData.length != A.getNumCols() )
throw new IllegalArgumentException("Unexpected sample length");
DenseMatrix64F mean = DenseMatrix64F.wrap(A.getNumCols(),1,this.mean);

DenseMatrix64F s = new DenseMatrix64F(A.getNumCols(),1,true,sampleData);
DenseMatrix64F r = new DenseMatrix64F(numComponents,1);

CommonOps.subtract(s, mean, s);

CommonOps.mult(V_t,s,r);

return r.data;
}

/**
* Converts a vector from eigen space into sample space.
*
* @param eigenData Eigen space data.
* @return Sample space projection.
*/
public double[] eigenToSampleSpace( double[] eigenData ) {
if( eigenData.length != numComponents )
throw new IllegalArgumentException("Unexpected sample length");

DenseMatrix64F s = new DenseMatrix64F(A.getNumCols(),1);
DenseMatrix64F r = DenseMatrix64F.wrap(numComponents,1,eigenData);

CommonOps.multTransA(V_t,r,s);

DenseMatrix64F mean = DenseMatrix64F.wrap(A.getNumCols(),1,this.mean);
CommonOps.add(s,mean,s);

return s.data;
}

/**
* 
* The membership error for a sample. If the error is less than a threshold then
* it can be considered a member. The threshold's value depends on the data set.
* 
* 
* The error is computed by projecting the sample into eigenspace then projecting
* it back into sample space and
* 
*
* @param sampleA The sample whose membership status is being considered.
* @return Its membership error.
*/
public double errorMembership( double[] sampleA ) {
double[] eig = sampleToEigenSpace(sampleA);
double[] reproj = eigenToSampleSpace(eig);

double total = 0;
for( int i = 0; i < reproj.length; i++ ) {
double d = sampleA[i] - reproj[i];
total += d*d;
}

return Math.sqrt(total);
}

/**
* Computes the dot product of each basis vector against the sample. Can be used as a measure
* for membership in the training sample set. High values correspond to a better fit.
*
* @param sample Sample of original data.
* @return Higher value indicates it is more likely to be a member of input dataset.
*/
public double response( double[] sample ) {
if( sample.length != A.numCols )
throw new IllegalArgumentException("Expected input vector to be in sample space");

DenseMatrix64F dots = new DenseMatrix64F(numComponents,1);
DenseMatrix64F s = DenseMatrix64F.wrap(A.numCols,1,sample);

CommonOps.mult(V_t,s,dots);

return NormOps.normF(dots);
}
}
</syntaxhighlight>

@@ Line 1: / Line 1: @@
 Principal Component Analysis (PCA) is a popular and simple to implement classification technique, often used in face recognition.  The following is an example of how to implement it in EJML using the procedural interface.  It is assumed that the reader is already familiar with PCA.
-Example on GitHub:
+External Resources
-* [https://github.com/lessthanoptimal/ejml/blob/v0.27/examples/src/org/ejml/example/PrincipalComponentAnalysis.java PrincipalComponentAnalysis]
+* [https://github.com/lessthanoptimal/ejml/blob/v0.27/examples/src/org/ejml/example/PrincipalComponentAnalysis.java PrincipalComponentAnalysis.java source code]
+* [http://en.wikipedia.org/wiki/Principal_component_analysis General PCA information on Wikipedia]
-For additional information on PCA:
+* <disqus>Discuss this example</disqus>
 = Sample Code =

@@ Line 2: / Line 2: @@
 External Resources
-* [https://github.com/lessthanoptimal/ejml/blob/v0.27/examples/src/org/ejml/example/PrincipalComponentAnalysis.java PrincipalComponentAnalysis.java source code]
+* [https://github.com/lessthanoptimal/ejml/blob/v0.31/examples/src/org/ejml/example/PrincipalComponentAnalysis.java PrincipalComponentAnalysis.java source code]
 * [http://en.wikipedia.org/wiki/Principal_component_analysis General PCA information on Wikipedia]
 * <disqus>Discuss this example</disqus>
@@ Line 40: / Line 40: @@
      // principal component subspace is stored in the rows
-     private DenseMatrix64F V_t;
+     private DMatrixRMaj V_t;
      // how many principal components are used
@@ Line 46: / Line 46: @@
      // where the data is stored
-     private DenseMatrix64F A = new DenseMatrix64F(1,1);
+     private DMatrixRMaj A = new DMatrixRMaj(1,1);
      private int sampleIndex;
@@ Line 120: / Line 120: @@
          // Compute SVD and save time by not computing U
-         SingularValueDecomposition<DenseMatrix64F> svd =
+         SingularValueDecomposition<DMatrixRMaj> svd =
-                 DecompositionFactory.svd(A.numRows, A.numCols, false, true, false);
+                 DecompositionFactory_DDRM.svd(A.numRows, A.numCols, false, true, false);
          if( !svd.decompose(A) )
              throw new RuntimeException("SVD failed");
          V_t = svd.getV(null,true);
-         DenseMatrix64F W = svd.getW(null);
+         DMatrixRMaj W = svd.getW(null);
          // Singular values are in an arbitrary order initially
-         SingularOps.descendingOrder(null,false,W,V_t,true);
+         SingularOps_DDRM.descendingOrder(null,false,W,V_t,true);
          // strip off unneeded components and find the basis
@@ Line 145: / Line 145: @@
              throw new IllegalArgumentException("Invalid component");
-         DenseMatrix64F v = new DenseMatrix64F(1,A.numCols);
+         DMatrixRMaj v = new DMatrixRMaj(1,A.numCols);
-         CommonOps.extract(V_t,which,which+1,0,A.numCols,v,0,0);
+         CommonOps_DDRM.extract(V_t,which,which+1,0,A.numCols,v,0,0);
          return v.data;
@@ Line 160: / Line 160: @@
          if( sampleData.length != A.getNumCols() )
              throw new IllegalArgumentException("Unexpected sample length");
-         DenseMatrix64F mean = DenseMatrix64F.wrap(A.getNumCols(),1,this.mean);
+         DMatrixRMaj mean = DMatrixRMaj.wrap(A.getNumCols(),1,this.mean);
-         DenseMatrix64F s = new DenseMatrix64F(A.getNumCols(),1,true,sampleData);
+         DMatrixRMaj s = new DMatrixRMaj(A.getNumCols(),1,true,sampleData);
-         DenseMatrix64F r = new DenseMatrix64F(numComponents,1);
+         DMatrixRMaj r = new DMatrixRMaj(numComponents,1);
-         CommonOps.subtract(s, mean, s);
+         CommonOps_DDRM.subtract(s, mean, s);
-         CommonOps.mult(V_t,s,r);
+         CommonOps_DDRM.mult(V_t,s,r);
          return r.data;
@@ Line 182: / Line 182: @@
              throw new IllegalArgumentException("Unexpected sample length");
-         DenseMatrix64F s = new DenseMatrix64F(A.getNumCols(),1);
+         DMatrixRMaj s = new DMatrixRMaj(A.getNumCols(),1);
-         DenseMatrix64F r = DenseMatrix64F.wrap(numComponents,1,eigenData);
+         DMatrixRMaj r = DMatrixRMaj.wrap(numComponents,1,eigenData);
-         CommonOps.multTransA(V_t,r,s);
+         CommonOps_DDRM.multTransA(V_t,r,s);
-         DenseMatrix64F mean = DenseMatrix64F.wrap(A.getNumCols(),1,this.mean);
+         DMatrixRMaj mean = DMatrixRMaj.wrap(A.getNumCols(),1,this.mean);
-         CommonOps.add(s,mean,s);
+         CommonOps_DDRM.add(s,mean,s);
          return s.data;
@@ Line 232: / Line 232: @@
              throw new IllegalArgumentException("Expected input vector to be in sample space");
-         DenseMatrix64F dots = new DenseMatrix64F(numComponents,1);
+         DMatrixRMaj dots = new DMatrixRMaj(numComponents,1);
-         DenseMatrix64F s = DenseMatrix64F.wrap(A.numCols,1,sample);
+         DMatrixRMaj s = DMatrixRMaj.wrap(A.numCols,1,sample);
-         CommonOps.mult(V_t,s,dots);
+         CommonOps_DDRM.mult(V_t,s,dots);
-         return NormOps.normF(dots);
+         return NormOps_DDRM.normF(dots);
+     }
+}
 </syntaxhighlight>