diff --git a/examples/correlation_matrices/sampler.cpp b/examples/correlation_matrices/sampler.cpp index 9a2ac08d2..f18f90a1b 100644 --- a/examples/correlation_matrices/sampler.cpp +++ b/examples/correlation_matrices/sampler.cpp @@ -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 MT; @@ -89,6 +90,28 @@ void write_to_file(std::string filename, std::vector 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 tol) + { + return false; + } + } + + //check if the matrix is positive semidefinite + using NT = double; + using MatrixType = Eigen::Matrix; + EigenvaluesProblems> solver; + + if(solver.isPositiveSemidefinite(matrix)) + { + return true; + } + return false; +} + template void correlation_matrix_uniform_sampling(const unsigned int n, const unsigned int num_points, std::string walkname){ @@ -106,7 +129,7 @@ void correlation_matrix_uniform_sampling(const unsigned int n, const unsigned in time = std::chrono::duration_cast(end - start).count(); std::cout << "Elapsed time : " << time << " (ms)" << std::endl; - write_to_file(walkname + "_matrices.txt", randPoints); + write_to_file(walkname + "_matrices" + std::to_string(n) + ".txt", randPoints); } template @@ -126,7 +149,15 @@ void correlation_matrix_uniform_sampling_MT(const unsigned int n, const unsigned time = std::chrono::duration_cast(end - start).count(); std::cout << "Elapsed time : " << time << " (ms)" << std::endl; - write_to_file(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(walkname + "_matrices_MT" + std::to_string(n) + ".txt", randPoints); } int main(int argc, char const *argv[]){ @@ -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 dimensions = {3, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100}; - old_uniform_sampling(n, num_points); + for(unsigned int n : dimensions){ - correlation_matrix_uniform_sampling(n, num_points, "BallWalk"); + old_uniform_sampling(n, num_points); - correlation_matrix_uniform_sampling(n, num_points, "RDHRWalk"); + correlation_matrix_uniform_sampling(n, num_points, "BallWalk"); - correlation_matrix_uniform_sampling(n, num_points, "BilliardWalk"); + correlation_matrix_uniform_sampling(n, num_points, "RDHRWalk"); - correlation_matrix_uniform_sampling(n, num_points, "AcceleratedBilliardWalk"); + correlation_matrix_uniform_sampling(n, num_points, "BilliardWalk"); - correlation_matrix_uniform_sampling_MT(n, num_points, "BallWalk"); + correlation_matrix_uniform_sampling(n, num_points, "AcceleratedBilliardWalk"); - correlation_matrix_uniform_sampling_MT(n, num_points, "RDHRWalk"); + correlation_matrix_uniform_sampling_MT(n, num_points, "BallWalk"); - correlation_matrix_uniform_sampling_MT(n, num_points, "BilliardWalk"); + correlation_matrix_uniform_sampling_MT(n, num_points, "RDHRWalk"); - correlation_matrix_uniform_sampling_MT(n, num_points, "AcceleratedBilliardWalk"); + correlation_matrix_uniform_sampling_MT(n, num_points, "BilliardWalk"); + + correlation_matrix_uniform_sampling_MT(n, num_points, "AcceleratedBilliardWalk"); + } return 0; -} \ No newline at end of file +} diff --git a/include/matrix_operations/EigenvaluesProblems.h b/include/matrix_operations/EigenvaluesProblems.h index 4012177f3..178d91666 100644 --- a/include/matrix_operations/EigenvaluesProblems.h +++ b/include/matrix_operations/EigenvaluesProblems.h @@ -167,9 +167,12 @@ class EigenvaluesProblems, E Spectra::DenseSymMatProd op(B); Spectra::DenseCholesky 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, Spectra::DenseCholesky, Spectra::GEIGS_CHOLESKY> - geigs(&op, &Bop, 1, 15 < matrixDim ? 15 : matrixDim); + geigs(&op, &Bop, 1, ncv); // Initialize and compute geigs.init(); @@ -324,9 +327,12 @@ class EigenvaluesProblems, E Spectra::DenseSymMatProd op(B); Spectra::DenseCholesky 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, Spectra::DenseCholesky, Spectra::GEIGS_CHOLESKY> - geigs(&op, &Bop, 1, 15 < matrixDim ? 15 : matrixDim); + geigs(&op, &Bop, 1, ncv); // Initialize and compute geigs.init(); @@ -439,10 +445,39 @@ class EigenvaluesProblems, 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 op(B); + Spectra::DenseCholesky Bop(A); + + //construct generalized eigen solver object, requesting the smallest eigenvalue + int ncv = std::min(std::max(10, matrixDim/20), matrixDim); + Spectra::SymGEigsSolver, Spectra::DenseCholesky, 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 ges(B,A); lambdaMinPositive = 1/ges.eigenvalues().reverse()[0]; eigvec = ges.eigenvectors().reverse().col(0).reverse(); +#endif return lambdaMinPositive; } };