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;
}
 
}  // namespace CASM
 
#endif