From f98178749ffcaf6aa42b3aba0a88b3f484942caf Mon Sep 17 00:00:00 2001 From: vboxuser Date: Fri, 31 May 2024 21:58:05 +0530 Subject: [PATCH 1/6] add spectra routine for minPosLinearEigenvalue_EigenSymSolver --- .../matrix_operations/EigenvaluesProblems.h | 30 ++++++++++++++++++- 1 file changed, 29 insertions(+), 1 deletion(-) diff --git a/include/matrix_operations/EigenvaluesProblems.h b/include/matrix_operations/EigenvaluesProblems.h index 4012177f3..62546d06b 100644 --- a/include/matrix_operations/EigenvaluesProblems.h +++ b/include/matrix_operations/EigenvaluesProblems.h @@ -324,7 +324,7 @@ 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 Spectra::SymGEigsSolver, Spectra::DenseCholesky, Spectra::GEIGS_CHOLESKY> geigs(&op, &Bop, 1, 15 < matrixDim ? 15 : matrixDim); @@ -439,10 +439,38 @@ 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(A); + Spectra::DenseCholesky Bop(B); + + //construct generalized eigen solver object, requesting the smallest eigenvalue + Spectra::SymGEigsSolver, Spectra::DenseCholesky, Spectra::GEIGS_CHOLESKY> + geigs(&op, &Bop, 1, 15 < matrixDim ? 15 : matrixDim); + + //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 = 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; } }; From 13ac65d91a624ccef0a59df2ffcea0c5793cf99a Mon Sep 17 00:00:00 2001 From: vboxuser Date: Tue, 4 Jun 2024 17:37:50 +0530 Subject: [PATCH 2/6] add new function to check for correlational matrices --- examples/correlation_matrices/sampler.cpp | 20 ++++++++++++++++++- .../matrix_operations/EigenvaluesProblems.h | 17 +++++++++------- .../random_walks/uniform_billiard_walk.hpp | 2 +- 3 files changed, 30 insertions(+), 9 deletions(-) diff --git a/examples/correlation_matrices/sampler.cpp b/examples/correlation_matrices/sampler.cpp index 9a2ac08d2..c903397eb 100644 --- a/examples/correlation_matrices/sampler.cpp +++ b/examples/correlation_matrices/sampler.cpp @@ -89,6 +89,24 @@ void write_to_file(std::string filename, std::vector const& randPoint std::cout.rdbuf(coutbuf); } +bool is_correlation_matrix(const MT& matrix){ + + //check if all the diagonal elements are ones + for(int i=0 ; i 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() > 0; +} + template void correlation_matrix_uniform_sampling(const unsigned int n, const unsigned int num_points, std::string walkname){ @@ -167,4 +185,4 @@ int main(int argc, char const *argv[]){ 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 62546d06b..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(); @@ -325,8 +326,9 @@ class EigenvaluesProblems, E Spectra::DenseCholesky Bop(-A); // 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(); @@ -443,12 +445,13 @@ class EigenvaluesProblems, E int matrixDim = A.rows(); NT lambdaMinPositive; - Spectra::DenseSymMatProd op(A); - Spectra::DenseCholesky Bop(B); + Spectra::DenseSymMatProd op(B); + Spectra::DenseCholesky Bop(A); //construct generalized eigen solver object, requesting the smallest eigenvalue - Spectra::SymGEigsSolver, Spectra::DenseCholesky, Spectra::GEIGS_CHOLESKY> - geigs(&op, &Bop, 1, 15 < matrixDim ? 15 : matrixDim); + 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(); @@ -462,7 +465,7 @@ class EigenvaluesProblems, E eigvec = geigs.eigenvectors().col(0); } - lambdaMinPositive = evalues(0); + lambdaMinPositive = NT(1)/evalues(0); #elif NT lambdaMinPositive = NT(0); diff --git a/include/random_walks/uniform_billiard_walk.hpp b/include/random_walks/uniform_billiard_walk.hpp index dcbc96b10..dd9528675 100644 --- a/include/random_walks/uniform_billiard_walk.hpp +++ b/include/random_walks/uniform_billiard_walk.hpp @@ -12,7 +12,7 @@ #define RANDOM_WALKS_UNIFORM_BILLIARD_WALK_HPP #include - +#include #include "convex_bodies/ball.h" #include "convex_bodies/ballintersectconvex.h" #include "convex_bodies/hpolytope.h" From 02fec522c258b73b83a469aa4da9a59bfcc4704a Mon Sep 17 00:00:00 2001 From: vboxuser Date: Wed, 5 Jun 2024 18:02:02 +0530 Subject: [PATCH 3/6] check for correlation matrices after sampling function --- examples/correlation_matrices/sampler.cpp | 21 ++++++++++++++++++--- 1 file changed, 18 insertions(+), 3 deletions(-) diff --git a/examples/correlation_matrices/sampler.cpp b/examples/correlation_matrices/sampler.cpp index c903397eb..1954b4c73 100644 --- a/examples/correlation_matrices/sampler.cpp +++ b/examples/correlation_matrices/sampler.cpp @@ -90,10 +90,12 @@ void write_to_file(std::string filename, std::vector const& randPoint } 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; } } @@ -104,7 +106,7 @@ bool is_correlation_matrix(const MT& matrix){ if(eigen_solver.info() != Eigen::Success) return false; //the matrix is positive definite if all eigenvalues are positive - return eigen_solver.eigenvalues().minCoeff() > 0; + return eigen_solver.eigenvalues().minCoeff() > tol; } template @@ -143,8 +145,21 @@ 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.txt", randPoints); + + //write_to_file(walkname + "_matrices_MT.txt", randPoints); } int main(int argc, char const *argv[]){ From 8501164d5c65084f2ddcdd49843ecdc590792b59 Mon Sep 17 00:00:00 2001 From: vboxuser Date: Thu, 6 Jun 2024 08:50:59 +0530 Subject: [PATCH 4/6] update tol in epsilon equality check --- examples/correlation_matrices/sampler.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/correlation_matrices/sampler.cpp b/examples/correlation_matrices/sampler.cpp index 1954b4c73..86882bf1d 100644 --- a/examples/correlation_matrices/sampler.cpp +++ b/examples/correlation_matrices/sampler.cpp @@ -106,7 +106,7 @@ bool is_correlation_matrix(const MT& 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; + return eigen_solver.eigenvalues().minCoeff() > -tol; } template From f2e5079bfee77cb6be061ce1a172af4a5b9d3cec Mon Sep 17 00:00:00 2001 From: vboxuser Date: Thu, 6 Jun 2024 09:42:19 +0530 Subject: [PATCH 5/6] check for matrices with greater dimensions --- examples/correlation_matrices/sampler.cpp | 31 +++++++++++++---------- 1 file changed, 17 insertions(+), 14 deletions(-) diff --git a/examples/correlation_matrices/sampler.cpp b/examples/correlation_matrices/sampler.cpp index 86882bf1d..b3bfa394c 100644 --- a/examples/correlation_matrices/sampler.cpp +++ b/examples/correlation_matrices/sampler.cpp @@ -126,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 @@ -157,9 +157,7 @@ void correlation_matrix_uniform_sampling_MT(const unsigned int n, const unsigned std::cout << "Number of valid points = " << valid_points << std::endl; - write_to_file(walkname + "_matrices_MT.txt", randPoints); - - //write_to_file(walkname + "_matrices_MT.txt", randPoints); + write_to_file(walkname + "_matrices_MT" + std::to_string(n) + ".txt", randPoints); } int main(int argc, char const *argv[]){ @@ -179,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; } From 9b7f114fcfe224c687d2f55c2d4903c414be477d Mon Sep 17 00:00:00 2001 From: atrayees <121290945+atrayees@users.noreply.github.com> Date: Thu, 6 Jun 2024 22:38:14 +0530 Subject: [PATCH 6/6] remove header header not required. --- include/random_walks/uniform_billiard_walk.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/include/random_walks/uniform_billiard_walk.hpp b/include/random_walks/uniform_billiard_walk.hpp index dd9528675..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 #include "convex_bodies/ball.h" #include "convex_bodies/ballintersectconvex.h" #include "convex_bodies/hpolytope.h"