pca.hpp

/*
 * eos - A 3D Morphable Model fitting library written in modern C++11/14.
 *
 * File: include/eos/pca/pca.hpp
 *
 * Copyright 2017 Patrik Huber
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 * http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
#pragma once
 
#ifndef PCA_HPP_
#define PCA_HPP_
 
#include "eos/morphablemodel/PcaModel.hpp"
 
#include "Eigen/Core"
#include "Eigen/Eigenvalues"
 
#include <array>
#include <cassert>
#include <utility>
#include <vector>
 
namespace eos {
namespace pca {
 
enum class Covariance {
    AtA, 
    AAt  
};
 
inline std::pair<Eigen::MatrixXf, Eigen::VectorXf> pca(const Eigen::Ref<const Eigen::MatrixXf> data,
                                                       Covariance covariance_type = Covariance::AtA)
{
    using Eigen::MatrixXf;
    using Eigen::VectorXf;
 
    MatrixXf cov;
    if (covariance_type == Covariance::AtA)
    {
        cov = data.adjoint() * data;
    } else if (covariance_type == Covariance::AAt)
    {
        cov = data * data.adjoint();
    }
 
    // The covariance is 1/(n-1) * AtA (or AAt), so divide by (num_samples - 1):
    cov /= (data.rows() - 1);
 
    const Eigen::SelfAdjointEigenSolver<MatrixXf> eig(cov);
 
    const auto num_eigenvectors_to_keep = data.rows() - 1;
 
    // Select the eigenvectors and eigenvalues that we want to keep, reverse them (from most significant to
    // least): For 'n' data points, we get at most 'n - 1' non-zero eigenvalues.
    VectorXf eigenvalues = eig.eigenvalues()
                               .bottomRows(num_eigenvectors_to_keep)
                               .reverse(); // eigenvalues() returns a column-vec
    MatrixXf eigenvectors = eig.eigenvectors().rightCols(num_eigenvectors_to_keep).rowwise().reverse();
 
    if (covariance_type == Covariance::AAt)
    {
        // Bring the AA^t variant in the right form by multiplying with A^t and 1/sqrt(eval):
        // (see e.g. https://math.stackexchange.com/questions/787822/how-do-covariance-matrix-c-aat-justify-the-actual-ata-in-pca)
        // (note the signs might be different from the AtA solution but that's not a problem as the sign of
        // eigenvectors are arbitrary anyway)
        eigenvectors = data.adjoint() * eigenvectors;
 
        // Multiply each eigenvector (column) with one over the square root of its respective eigenvalue
        // (1/sqrt(eigenvalue(i))): (this is a neat short-hand notation, see
        // https://stackoverflow.com/a/42945996/1345959).
        const VectorXf one_over_sqrt_eigenvalues = eigenvalues.array().rsqrt();
        eigenvectors *= one_over_sqrt_eigenvalues.asDiagonal();
 
        // Compensate for the covariance division by (n - 1) above:
        eigenvectors /= std::sqrt(data.rows() - 1);
    }
 
    return { eigenvectors, eigenvalues };
};
 
inline std::pair<Eigen::MatrixXf, Eigen::VectorXf> pca(const Eigen::Ref<const Eigen::MatrixXf> data,
                                                       int num_eigenvectors_to_keep,
                                                       Covariance covariance_type = Covariance::AtA)
{
    using Eigen::MatrixXf;
    using Eigen::VectorXf;
 
    VectorXf eigenvalues;
    MatrixXf eigenvectors;
    std::tie(eigenvectors, eigenvalues) = pca(data, covariance_type);
 
    // Reduce the basis and eigenvalues, and return:
    assert(num_eigenvectors_to_keep <= eigenvectors.size());
    return { eigenvectors.leftCols(num_eigenvectors_to_keep), eigenvalues.topRows(num_eigenvectors_to_keep) };
};
 
inline std::pair<Eigen::MatrixXf, Eigen::VectorXf> pca(const Eigen::Ref<const Eigen::MatrixXf> data,
                                                       float variance_to_keep,
                                                       Covariance covariance_type = Covariance::AtA)
{
    using Eigen::MatrixXf;
    using Eigen::VectorXf;
 
    assert(variance_to_keep >= 0.0f && variance_to_keep <= 1.0f);
 
    VectorXf eigenvalues;
    MatrixXf eigenvectors;
    std::tie(eigenvectors, eigenvalues) = pca(data, covariance_type);
 
    // Figure out how many coeffs to keep:
    // variance_explained_by_first_comp = eigenval(1)/sum(eigenvalues)
    // variance_explained_by_second_comp = eigenval(2)/sum(eigenvalues), etc.
    auto num_eigenvectors_to_keep = eigenvalues.size(); // In the "worst" case we return all eigenvectors.
    const auto total_sum = eigenvalues.sum();
    float cum_sum = 0.0f;
    for (int i = 0; i < eigenvalues.size(); ++i)
    {
        cum_sum += eigenvalues(i);
        // If the current variation explained is larger or equal to the amount of variation that
        // the user requested to keep, we're done:
        if (cum_sum / total_sum >= variance_to_keep)
        {
            num_eigenvectors_to_keep = i + 1;
            break;
        }
    }
 
    // Reduce the basis and eigenvalues, and return:
    assert(num_eigenvectors_to_keep <= eigenvectors.size());
    return { eigenvectors.leftCols(num_eigenvectors_to_keep), eigenvalues.topRows(num_eigenvectors_to_keep) };
};
 
inline morphablemodel::PcaModel pca(const Eigen::Ref<const Eigen::MatrixXf> data,
                                    std::vector<std::array<int, 3>> triangle_list,
                                    Covariance covariance_type = Covariance::AtA)
{
    using Eigen::MatrixXf;
    using Eigen::VectorXf;
 
    // Compute the mean and mean-free data matrix:
    // Each row is one instance of data (e.g. a 3D scan)
    const VectorXf mean = data.colwise().mean();
    const MatrixXf meanfree_data = data.rowwise() - mean.transpose();
 
    // Carry out PCA and return the constructed model:
    VectorXf eigenvalues;
    MatrixXf eigenvectors;
    std::tie(eigenvectors, eigenvalues) = pca(meanfree_data, covariance_type);
 
    return morphablemodel::PcaModel(mean, eigenvectors, eigenvalues, triangle_list);
};
 
} /* namespace pca */
} /* namespace eos */
 
#endif /* PCA_HPP_ */