Skip to content

Commit

Permalink
Feature/spectra correlations (#306)
Browse files Browse the repository at this point in the history
  • Loading branch information
@atrayees
atrayees authored Jun 17, 2024
1 parent 3129f7b commit e484865
Show file tree
Hide file tree
Showing 2 changed files with 88 additions and 17 deletions.
62 changes: 49 additions & 13 deletions examples/correlation_matrices/sampler.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@
#include "convex_bodies/spectrahedra/spectrahedron.h"
#include "random_walks/random_walks.hpp"
#include "sampling/sample_correlation_matrices.hpp"
#include "matrix_operations/EigenvaluesProblems.h"

typedef double NT;
typedef Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic> MT;
Expand Down Expand Up @@ -89,6 +90,28 @@ void write_to_file(std::string filename, std::vector<PointType> const& randPoint
std::cout.rdbuf(coutbuf);
}

bool is_correlation_matrix(const MT& matrix, const double tol = 1e-8){
//check if all the diagonal elements are ones
for(int i=0 ; i<matrix.rows() ; i++)
{
if(std::abs(matrix(i, i)-1.0) > tol)
{
return false;
}
}

//check if the matrix is positive semidefinite
using NT = double;
using MatrixType = Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic>;
EigenvaluesProblems<NT, MatrixType, Eigen::Matrix<NT, Eigen::Dynamic, 1>> solver;

if(solver.isPositiveSemidefinite(matrix))
{
return true;
}
return false;
}

template<typename WalkType>
void correlation_matrix_uniform_sampling(const unsigned int n, const unsigned int num_points, std::string walkname){

Expand All @@ -106,7 +129,7 @@ void correlation_matrix_uniform_sampling(const unsigned int n, const unsigned in
time = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
std::cout << "Elapsed time : " << time << " (ms)" << std::endl;

write_to_file<Point>(walkname + "_matrices.txt", randPoints);
write_to_file<Point>(walkname + "_matrices" + std::to_string(n) + ".txt", randPoints);
}

template<typename WalkType>
Expand All @@ -126,7 +149,15 @@ void correlation_matrix_uniform_sampling_MT(const unsigned int n, const unsigned
time = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
std::cout << "Elapsed time : " << time << " (ms)" << std::endl;

write_to_file<PointMT>(walkname + "_matrices_MT.txt", randPoints);
int valid_points = 0;
for(const auto& points : randPoints){
if(is_correlation_matrix(points.mat)){
valid_points++;
}
}
std::cout << "Number of valid points = " << valid_points << std::endl;

write_to_file<PointMT>(walkname + "_matrices_MT" + std::to_string(n) + ".txt", randPoints);
}

int main(int argc, char const *argv[]){
Expand All @@ -146,25 +177,30 @@ int main(int argc, char const *argv[]){
printf("\n");
#endif

unsigned int n = 3, num_points = 5000;
unsigned int num_points = 5000;

std::vector<unsigned int> dimensions = {3, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100};

old_uniform_sampling<BilliardWalk>(n, num_points);
for(unsigned int n : dimensions){

correlation_matrix_uniform_sampling<BallWalk>(n, num_points, "BallWalk");
old_uniform_sampling<BilliardWalk>(n, num_points);

correlation_matrix_uniform_sampling<RDHRWalk>(n, num_points, "RDHRWalk");
correlation_matrix_uniform_sampling<BallWalk>(n, num_points, "BallWalk");

correlation_matrix_uniform_sampling<BilliardWalk>(n, num_points, "BilliardWalk");
correlation_matrix_uniform_sampling<RDHRWalk>(n, num_points, "RDHRWalk");

correlation_matrix_uniform_sampling<AcceleratedBilliardWalk>(n, num_points, "AcceleratedBilliardWalk");
correlation_matrix_uniform_sampling<BilliardWalk>(n, num_points, "BilliardWalk");

correlation_matrix_uniform_sampling_MT<BallWalk>(n, num_points, "BallWalk");
correlation_matrix_uniform_sampling<AcceleratedBilliardWalk>(n, num_points, "AcceleratedBilliardWalk");

correlation_matrix_uniform_sampling_MT<RDHRWalk>(n, num_points, "RDHRWalk");
correlation_matrix_uniform_sampling_MT<BallWalk>(n, num_points, "BallWalk");

correlation_matrix_uniform_sampling_MT<BilliardWalk>(n, num_points, "BilliardWalk");
correlation_matrix_uniform_sampling_MT<RDHRWalk>(n, num_points, "RDHRWalk");

correlation_matrix_uniform_sampling_MT<AcceleratedBilliardWalk>(n, num_points, "AcceleratedBilliardWalk");
correlation_matrix_uniform_sampling_MT<BilliardWalk>(n, num_points, "BilliardWalk");

correlation_matrix_uniform_sampling_MT<AcceleratedBilliardWalk>(n, num_points, "AcceleratedBilliardWalk");
}

return 0;
}
}
43 changes: 39 additions & 4 deletions include/matrix_operations/EigenvaluesProblems.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -167,9 +167,12 @@ class EigenvaluesProblems<NT, Eigen::Matrix<NT,Eigen::Dynamic,Eigen::Dynamic>, E
Spectra::DenseSymMatProd<NT> op(B);
Spectra::DenseCholesky<NT> Bop(-A);

// Construct generalized eigen solver object, requesting the largest three generalized eigenvalues
// Construct generalized eigen solver object, computing the minimum positive eigenvalue by computing the largest eigenvalue of the inverse Generalized Eigenvalue Problem
// An empirical value of ncv that gives a better performance
// TODO: tune this implementation by tuning the parameters like ncv
int ncv = std::min(std::max(10, matrixDim/20), matrixDim);
Spectra::SymGEigsSolver<NT, Spectra::LARGEST_ALGE, Spectra::DenseSymMatProd<NT>, Spectra::DenseCholesky<NT>, Spectra::GEIGS_CHOLESKY>
geigs(&op, &Bop, 1, 15 < matrixDim ? 15 : matrixDim);
geigs(&op, &Bop, 1, ncv);

// Initialize and compute
geigs.init();
Expand Down Expand Up @@ -324,9 +327,12 @@ class EigenvaluesProblems<NT, Eigen::Matrix<NT,Eigen::Dynamic,Eigen::Dynamic>, E
Spectra::DenseSymMatProd<NT> op(B);
Spectra::DenseCholesky<NT> Bop(-A);

// Construct generalized eigen solver object, requesting the largest three generalized eigenvalues
// Construct generalized eigen solver object, requesting the largest generalized eigenvalue
// an empirical value of ncv that gives a better performance
// TODO: tune this implementation by tuning the parameters like ncv
int ncv = std::min(std::max(10, matrixDim/20), matrixDim);
Spectra::SymGEigsSolver<NT, Spectra::LARGEST_ALGE, Spectra::DenseSymMatProd<NT>, Spectra::DenseCholesky<NT>, Spectra::GEIGS_CHOLESKY>
geigs(&op, &Bop, 1, 15 < matrixDim ? 15 : matrixDim);
geigs(&op, &Bop, 1, ncv);

// Initialize and compute
geigs.init();
Expand Down Expand Up @@ -439,10 +445,39 @@ class EigenvaluesProblems<NT, Eigen::Matrix<NT,Eigen::Dynamic,Eigen::Dynamic>, E
/// \param[in] B: symmetric matrix
/// \return The minimum positive eigenvalue and the corresponding eigenvector
NT minPosLinearEigenvalue_EigenSymSolver(MT const & A, MT const & B, VT &eigvec) const {

#if defined(SPECTRA_EIGENVALUES_SOLVER)
int matrixDim = A.rows();
NT lambdaMinPositive;

Spectra::DenseSymMatProd<NT> op(B);
Spectra::DenseCholesky<NT> Bop(A);

//construct generalized eigen solver object, requesting the smallest eigenvalue
int ncv = std::min(std::max(10, matrixDim/20), matrixDim);
Spectra::SymGEigsSolver<NT, Spectra::LARGEST_ALGE, Spectra::DenseSymMatProd<NT>, Spectra::DenseCholesky<NT>, Spectra::GEIGS_CHOLESKY>
geigs(&op, &Bop, 1, ncv);

//initialize and compute
geigs.init();
int nconv = geigs.compute();

//retrieve results
VT evalues;

if(geigs.info() == Spectra::SUCCESSFUL){
evalues = geigs.eigenvalues();
eigvec = geigs.eigenvectors().col(0);
}

lambdaMinPositive = NT(1)/evalues(0);

#elif
NT lambdaMinPositive = NT(0);
Eigen::GeneralizedSelfAdjointEigenSolver<MT> ges(B,A);
lambdaMinPositive = 1/ges.eigenvalues().reverse()[0];
eigvec = ges.eigenvectors().reverse().col(0).reverse();
#endif
return lambdaMinPositive;
}
};
Expand Down

0 comments on commit e484865

Please sign in to comment.