From 1c110354a828bc49e3cef8060118cc9800bbc951 Mon Sep 17 00:00:00 2001 From: Vissarion Fisikopoulos Date: Fri, 7 Jun 2024 12:26:03 +0300 Subject: [PATCH] Revert "add spectra routine for minPosLinearEigenvalue_EigenSymSolver (#303)" This reverts commit 9a6676a1ce7d4e6dde28fd290e11b575d88ff837. --- examples/correlation_matrices/sampler.cpp | 64 ++++--------------- .../matrix_operations/EigenvaluesProblems.h | 37 +---------- .../random_walks/uniform_billiard_walk.hpp | 1 + 3 files changed, 18 insertions(+), 84 deletions(-) diff --git a/examples/correlation_matrices/sampler.cpp b/examples/correlation_matrices/sampler.cpp index b3bfa394c..9a2ac08d2 100644 --- a/examples/correlation_matrices/sampler.cpp +++ b/examples/correlation_matrices/sampler.cpp @@ -89,26 +89,6 @@ 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){ @@ -126,7 +106,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" + std::to_string(n) + ".txt", randPoints); + write_to_file(walkname + "_matrices.txt", randPoints); } template @@ -145,19 +125,8 @@ 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; - - 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); + + write_to_file(walkname + "_matrices_MT.txt", randPoints); } int main(int argc, char const *argv[]){ @@ -177,30 +146,25 @@ int main(int argc, char const *argv[]){ printf("\n"); #endif - unsigned int num_points = 5000; - - std::vector dimensions = {3, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100}; - - for(unsigned int n : dimensions){ + unsigned int n = 3, num_points = 5000; - 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 1beea5d2c..4012177f3 100644 --- a/include/matrix_operations/EigenvaluesProblems.h +++ b/include/matrix_operations/EigenvaluesProblems.h @@ -168,9 +168,8 @@ 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, ncv); + geigs(&op, &Bop, 1, 15 < matrixDim ? 15 : matrixDim); // Initialize and compute geigs.init(); @@ -325,10 +324,9 @@ class EigenvaluesProblems, E Spectra::DenseSymMatProd op(B); Spectra::DenseCholesky Bop(-A); - // Construct generalized eigen solver object, requesting the largest generalized eigenvalue - int ncv = std::min(std::max(10, matrixDim/20), matrixDim); + // Construct generalized eigen solver object, requesting the largest three generalized eigenvalues Spectra::SymGEigsSolver, Spectra::DenseCholesky, Spectra::GEIGS_CHOLESKY> - geigs(&op, &Bop, 1, ncv); + geigs(&op, &Bop, 1, 15 < matrixDim ? 15 : matrixDim); // Initialize and compute geigs.init(); @@ -441,39 +439,10 @@ 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 14729d27d..dcbc96b10 100644 --- a/include/random_walks/uniform_billiard_walk.hpp +++ b/include/random_walks/uniform_billiard_walk.hpp @@ -12,6 +12,7 @@ #define RANDOM_WALKS_UNIFORM_BILLIARD_WALK_HPP #include + #include "convex_bodies/ball.h" #include "convex_bodies/ballintersectconvex.h" #include "convex_bodies/hpolytope.h"