PCA.hh

Go to the documentation of this file.

#ifndef CASM_PCA
#define CASM_PCA

#include "casm/external/Eigen/Dense"

namespace CASM {

  class PCA {

  public:

    typedef Eigen::MatrixXd::Index size_type;


    template<typename Derived>
    PCA(const Eigen::MatrixBase<Derived> &M, double singular_value_tol = 1e-14) {

      // set the mean to zero
      Eigen::MatrixXd mat_mean_zero = M;
      mat_mean_zero.colwise() -= (M.rowwise().mean());

      // principle component analysis to rotate to input space that to its range
      Eigen::JacobiSVD<Eigen::MatrixXd> svd(mat_mean_zero, Eigen::ComputeFullU);
      int rank = 0;
      for(Index i = 0; i < svd.singularValues().size(); i++) {
        if(std::abs(svd.singularValues()[i]) <= singular_value_tol) {
          break;
        }
        else {
          rank++;
        }
      }

      m_reduce = svd.matrixU().block(0, 0, M.rows(), rank).transpose();

    }

    size_type dim() const {
      return m_reduce.cols();
    }

    size_type rank() const {
      return m_reduce.rows();
    }

    Eigen::MatrixXd reduce() const {
      return m_reduce;
    }

    Eigen::MatrixXd expand() const {
      return m_reduce.transpose();
    }


  private:

    // transform full dimension [input...] column vector onto subspace [range(input...)]
    Eigen::MatrixXd m_reduce;

  };

  inline Eigen::MatrixXd pad(const Eigen::MatrixXd &M, int n) {
    Eigen::MatrixXd result = Eigen::MatrixXd::Zero(M.rows() + n, M.cols() + n);
    result << M, Eigen::MatrixXd::Zero(M.rows(), n),
           Eigen::MatrixXd::Zero(n, M.cols()), Eigen::MatrixXd::Identity(n, n);
    return result;
  }

}


#endif