pca.hpp

Go to the documentation of this file.

#ifndef IRLBA_PCA_HPP
#define IRLBA_PCA_HPP
 
#include <cmath>
#include <utility>
#include <stdexcept>
 
#include "compute.hpp"
#include "Matrix/simple.hpp"
#include "Matrix/centered.hpp"
#include "Matrix/scaled.hpp"
 
#include "Eigen/Dense"
 
namespace irlba {
 
template<class InputEigenMatrix_, class OutputEigenMatrix_, class EigenVector_>
std::pair<bool, int> pca(
    const InputEigenMatrix_& matrix,
    bool center,
    bool scale,
    Eigen::Index number,
    OutputEigenMatrix_& scores,
    OutputEigenMatrix_& rotation,
    EigenVector_& variances,
    const Options<EigenVector_>& options
) {
    const Eigen::Index nr = matrix.rows();
    const Eigen::Index nc = matrix.cols();
 
    EigenVector_ center0;
    if (center) {
        if (nr < 1) {
            throw std::runtime_error("cannot center with no observations");    
        }
        center0.resize(nc);
    }
 
    EigenVector_ scale0;
    if (scale) {
        if (nr < 2) {
            throw std::runtime_error("cannot scale with fewer than two observations");    
        }
        scale0.resize(nc);
    }
 
    if (center || scale) {
        EigenVector_ buffer;
        for (Eigen::Index i = 0; i < nc; ++i) {
            buffer = matrix.col(i);
 
            typename EigenVector_::Scalar mean = 0;
            if (center) {
                mean = buffer.sum() / nr;
                center0[i] = mean;
            }
 
            if (scale) {
                typename EigenVector_::Scalar var = 0;
                for (const auto x : buffer) {
                    var += (x - mean) * (x - mean);
                }
 
                if (var) {
                    scale0[i] = std::sqrt(var/(nr - 1));
                } else {
                    scale0[i] = 1;
                }
            }
        }
    }
 
    // Reduce the size of the binary by forcing all compute() calls to use virtual dispatch,
    // rather than realizing four new instances of the same function with different Matrix subclasses.
    std::unique_ptr<Matrix<EigenVector_, OutputEigenMatrix_> > ptr;
    ptr.reset(new SimpleMatrix<EigenVector_, OutputEigenMatrix_, I<decltype(&matrix)> >(&matrix));
 
    if (center) {
        std::unique_ptr<Matrix<EigenVector_, OutputEigenMatrix_> > alt;
        alt.reset(new CenteredMatrix<EigenVector_, OutputEigenMatrix_, I<decltype(ptr)>, I<decltype(&center0)> >(std::move(ptr), &center0));
        ptr.swap(alt);
    }
 
    if (scale) {
        std::unique_ptr<Matrix<EigenVector_, OutputEigenMatrix_> > alt;
        alt.reset(new ScaledMatrix<EigenVector_, OutputEigenMatrix_, I<decltype(ptr)>, I<decltype(&scale0)> >(std::move(ptr), &scale0, true, true));
        ptr.swap(alt);
    }
 
    const auto stats = compute(*ptr, number, scores, rotation, variances, options);
 
    scores.array().rowwise() *= variances.adjoint().array();
    if (nc > 1) {
        const auto denom = nr - 1;
        for (auto& v : variances) {
            v = v * v / denom;
        }
    }
 
    return stats;
}
 
template<class EigenMatrix_, class EigenVector_>
struct PcaResults {
    EigenMatrix_ scores;
 
    EigenMatrix_ rotation;
 
    EigenVector_ variances;
 
    int iterations;
 
    bool converged;
};
 
template<class OutputEigenMatrix_ = Eigen::MatrixXd, class EigenVector_ = Eigen::VectorXd, class InputEigenMatrix_>
PcaResults<OutputEigenMatrix_, EigenVector_> pca(const InputEigenMatrix_& matrix, bool center, bool scale, Eigen::Index number, const Options<EigenVector_>& options) {
    PcaResults<OutputEigenMatrix_, EigenVector_> output;
    const auto stats = pca(matrix, center, scale, number, output.scores, output.rotation, output.variances, options);
    output.converged = stats.first;
    output.iterations = stats.second;
    return output;
}
 
}
 
#endif