diff --git a/examples/correlation_matrices/sampler.cpp b/examples/correlation_matrices/sampler.cpp index 9a2ac08d2..b3bfa394c 100644 --- a/examples/correlation_matrices/sampler.cpp +++ b/examples/correlation_matrices/sampler.cpp @@ -89,6 +89,26 @@ 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 definite + Eigen::SelfAdjointEigenSolver eigen_solver(matrix); + + if(eigen_solver.info() != Eigen::Success) return false; + + //the matrix is positive definite if all eigenvalues are positive + return eigen_solver.eigenvalues().minCoeff() > -tol; +} + template void correlation_matrix_uniform_sampling(const unsigned int n, const unsigned int num_points, std::string walkname){ @@ -106,7 +126,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 @@ -125,8 +145,19 @@ void correlation_matrix_uniform_sampling_MT(const unsigned int n, const unsigned end = std::chrono::steady_clock::now(); 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}; + + for(unsigned int n : dimensions){ - old_uniform_sampling(n, num_points); + old_uniform_sampling(n, num_points); - correlation_matrix_uniform_sampling(n, num_points, "BallWalk"); + correlation_matrix_uniform_sampling(n, num_points, "BallWalk"); - correlation_matrix_uniform_sampling(n, num_points, "RDHRWalk"); + correlation_matrix_uniform_sampling(n, num_points, "RDHRWalk"); - correlation_matrix_uniform_sampling(n, num_points, "BilliardWalk"); + correlation_matrix_uniform_sampling(n, num_points, "BilliardWalk"); - correlation_matrix_uniform_sampling(n, num_points, "AcceleratedBilliardWalk"); + correlation_matrix_uniform_sampling(n, num_points, "AcceleratedBilliardWalk"); - correlation_matrix_uniform_sampling_MT(n, num_points, "BallWalk"); + correlation_matrix_uniform_sampling_MT(n, num_points, "BallWalk"); - correlation_matrix_uniform_sampling_MT(n, num_points, "RDHRWalk"); + correlation_matrix_uniform_sampling_MT(n, num_points, "RDHRWalk"); - correlation_matrix_uniform_sampling_MT(n, num_points, "BilliardWalk"); + correlation_matrix_uniform_sampling_MT(n, num_points, "BilliardWalk"); - correlation_matrix_uniform_sampling_MT(n, num_points, "AcceleratedBilliardWalk"); + 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..1beea5d2c 100644 --- a/include/matrix_operations/EigenvaluesProblems.h +++ b/include/matrix_operations/EigenvaluesProblems.h @@ -168,8 +168,9 @@ class EigenvaluesProblems, E Spectra::DenseCholesky Bop(-A); // Construct generalized eigen solver object, requesting the largest three generalized eigenvalues + 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 +325,10 @@ 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 + 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 +441,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; } }; diff --git a/include/random_walks/uniform_billiard_walk.hpp b/include/random_walks/uniform_billiard_walk.hpp index dcbc96b10..14729d27d 100644 --- a/include/random_walks/uniform_billiard_walk.hpp +++ b/include/random_walks/uniform_billiard_walk.hpp @@ -12,7 +12,6 @@ #define RANDOM_WALKS_UNIFORM_BILLIARD_WALK_HPP #include - #include "convex_bodies/ball.h" #include "convex_bodies/ballintersectconvex.h" #include "convex_bodies/hpolytope.h"