rtcPCA.h

Go to the documentation of this file.

/*
 * Copyright (C) 2008
 * Robert Bosch LLC
 * Research and Technology Center North America
 * Palo Alto, California
 *
 * All rights reserved.
 *
 *------------------------------------------------------------------------------
 * project ....: PUMA: Probablistic Unsupervised Model Acquisition
 * file .......: pca.h
 * authors ....: Benjamin Pitzer
 * organization: Robert Bosch LLC
 * creation ...: 08/16/2006
 * modified ...: $Date: 2009-01-10 22:38:24 -0800 (Sat, 10 Jan 2009) $
 * changed by .: $Author: benjaminpitzer $
 * revision ...: $Revision: 695 $
 */
#ifndef PCA_H
#define PCA_H

#include "rtc/rtcMath.h"
#include "rtc/rtcSMat.h"
#include "rtc/rtcVec.h"
#include "rtc/rtcVarMat.h"
#include "rtc/rtcVarVec.h"
#include "rtc/rtcEigenSystem.h"

namespace rtc {

  template <class T, int N>
  class PrincipalComponentAnalysis {
  public:
    // data vector type
    typedef std::vector< Vec<T,N> > VecData;

    PrincipalComponentAnalysis();
    ~PrincipalComponentAnalysis();

    SMat<T,N>& computeTransformMatrix(const VecData& data);
    SMat<T,N>& computeTransformMatrix(const VarMat<T>& data);

    Vec<T,N>& transform(const Vec<T,N>& src, Vec<T,N>& dest) const;
    Vec<T,N>& transform(Vec<T,N>& srcdest) const;
    VarMat<T>& transform(VarMat<T>& data) const;
    VecData& transform(VecData& data) const;
  protected:
    // calculates mean of all rows
    SMat<T,N>& covarianceMatrixOfRows(const VecData& src, SMat<T,N>& covariance) const;
    // calculates mean of all rows
    Vec<T,N>& meanOfRows(const VecData& src, Vec<T,N>& mean) const;

    // Data
    bool transform_matrix_computed;

  public:
    Vec<T,N> offset;
    SMat<T,N> correlation_coefficient;
    SMat<T,N> eigen_vectors;
    SMat<T,N> transform_matrix;
    Vec<T,N> eigen_values;
  };


  template <class T, int N>
  inline PrincipalComponentAnalysis<T,N>::PrincipalComponentAnalysis() {
    transform_matrix_computed = false;
  }

  template <class T, int N>
  inline PrincipalComponentAnalysis<T,N>::~PrincipalComponentAnalysis() {

  }

  template <class T, int N>
  inline SMat<T,N>& PrincipalComponentAnalysis<T,N>::computeTransformMatrix(const VecData& data)
  {
    // determine the center of the data point distribution
    meanOfRows(data,offset);

    // correlation coefficient matrix
    covarianceMatrixOfRows(data,correlation_coefficient);

    // now compute eigenvectors of correlation coefficient matrix
    EigenSystem<T,N> eig(correlation_coefficient);
    eig.jacobi(eigen_values, eigen_vectors);

    // sort eigenvalues (and corresponding eigenvectors) into descending order
    eig.eigenSort(eigen_values, eigen_vectors);

    // the transformation matrix is simply the first columns of
    // the ordered eigenvectors
    transform_matrix.set(eigen_vectors);

    transform_matrix_computed = true;
    return transform_matrix;
  }

  template <class T, int N>
  inline SMat<T,N>& PrincipalComponentAnalysis<T,N>::computeTransformMatrix(const VarMat<T>& data)
  {
    if(data.columns()!=N)
      throw Exception("The input matrix has the wrong number of columns.");

    // determine the center of the data point distribution
    VarVec<T> var_offset = data.meanOfRows();
    for(int i=0;i<N;++i) offset(i)=var_offset(i);

    // correlation coefficient matrix
    VarSMat<T> var_correlation_coefficient = data.covarianceMatrixOfRows();
    for(int i=0;i<N;++i) for(int j=0;j<N;++j) correlation_coefficient(i,j)=var_correlation_coefficient(i,j);

    // now compute eigenvectors of correlation coefficient matrix
    EigenSystem<T,N> eig(correlation_coefficient);
    eig.jacobi(eigen_values, eigen_vectors);

    // sort eigenvalues (and corresponding eigenvectors) into descending order
    eig.eigenSort(eigen_values, eigen_vectors);

    // the transformation matrix is simply the first columns of
    // the ordered eigenvectors
    transform_matrix.set(eigen_vectors);

    transform_matrix_computed = true;
    return transform_matrix;
  }

  /*
   * Transforms a single vector according to a previously computed
   * transformation matrix.
   */
  template <class T, int N>
  inline Vec<T,N>& PrincipalComponentAnalysis<T,N>::transform(const Vec<T,N>& src, Vec<T,N>& dest) const {
    // check if transform matrix was computed
    if(!transform_matrix_computed)
      throw Exception("The transform matrix has to be computed before calling this function.");
    // subtract offset
    Vec<T,N> tmp = src-offset;
    // apply transform
    dest = transform_matrix*tmp;
    return dest;
  }

  /*
   * Transforms a single vector according to a previously computed
   * transformation matrix.
   */
  template <class T, int N>
  inline Vec<T,N>& PrincipalComponentAnalysis<T,N>::transform(Vec<T,N>& srcdest) const {
    // check if transform matrix was computed
    if(!transform_matrix_computed)
      throw Exception("The transform matrix has to be computed before calling this function.");
    // subtract offset
    Vec<T,N> tmp = srcdest-offset;
    // apply transform
    srcdest = transform_matrix*tmp;
    return srcdest;
  }

  /*
   * Transforms the data vector according to a previously computed
   * transformation matrix.
   */
  template <class T, int N>
  inline VarMat<T>& PrincipalComponentAnalysis<T,N>::transform(VarMat<T>& data) const {
    // check if transform matrix was computed
    if(!transform_matrix_computed)
      throw Exception("The transform matrix has to be computed before calling this function.");

    for (unsigned int i=0;i<data.rows();i++) {
      Vec<T,N> tmp = data.getRow(i);
      transform(tmp);
      data.setRow(i,tmp);
    }

    return data;
  }

  /*
   * Transforms the data vector according to a previously computed
   * transformation matrix.
   */
  template <class T, int N>
  inline std::vector< Vec<T,N> >& PrincipalComponentAnalysis<T,N>::transform(VecData& data) const {
    // check if transform matrix was computed
    if(!transform_matrix_computed)
      throw Exception("The transform matrix has to be computed before calling this function.");

    for (unsigned int i=0;i<data.size();i++)
      transform(data[i]);

    return data;
  }

  template <class T, int N>
  inline Vec<T,N>& PrincipalComponentAnalysis<T,N>::meanOfRows(const VecData& src, Vec<T,N>& mean) const {
    mean.set(T(0));
    int len = src.size();
    for (int j=0;j<N;j++)
      for (int i=0;i<len;i++) mean.x[j] += src[i][j];
    mean/=static_cast<T>(len);
    return mean;
  }

  // covariance of rows
  template <class T, int N>
  inline SMat<T,N>& PrincipalComponentAnalysis<T,N>::covarianceMatrixOfRows(const VecData& src, SMat<T,N>& covariance) const {

    int len = src.size();
    covariance.set(T(0));

    //Implementation for the sum of Matrix[X];
    if (len<2) { //just one row
      throw Exception("The input data has to have more than one row.");
    }

    // mean vector
    Vec<T,N> mu(T(0));

    /*
     * the loop gets out the result of the Matrix[X];
     * Matrix[X]=sum(Matrix[X(i)]);-- i from 0 to n-1--
     */
    for (int i=0; i<len; i++) {
      const Vec<T,N>& tmpV = src.at(i);
      mu.add(tmpV);
      SMat<T,N> tmp = tmpV.outer(tmpV);
      covariance.add(tmp); //add the Matrix[X(i)] to Matrix[sumMatrix];
    }
    //End of the Implementation for the sum of Matrix[X];

    //Implementation for the Matrix[meanOf];
    mu/=static_cast<T>(len);

    SMat<T,N> tmp = mu.outer(mu);
    tmp *= static_cast<T>(len);

    //get the result of the difference from Matrix[X] and Matrix[meanOf];
    covariance.subtract(tmp);
    covariance/=static_cast<T>(len-1);
    return covariance;
  }

}; // end namspace puma

#endif