Skip to content

Commit

Permalink
add spectra routine for minPosLinearEigenvalue_EigenSymSolver (#303)
Browse files Browse the repository at this point in the history 
* add spectra routine for minPosLinearEigenvalue_EigenSymSolver

* add new function to check for correlational matrices

* check for correlation matrices after sampling function

* update tol in epsilon equality check

* check for matrices with greater dimensions

* remove <iostream> header

header not required.

---------

Co-authored-by: vboxuser <[email protected]>
  • Loading branch information
@atrayees @kasryan
atrayees and kasryan authored Jun 6, 2024
1 parent 6a5a17e commit 9a6676a
Show file tree
Hide file tree
Showing 3 changed files with 84 additions and 18 deletions.
64 changes: 50 additions & 14 deletions examples/correlation_matrices/sampler.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -89,6 +89,26 @@ 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 definite
Eigen::SelfAdjointEigenSolver<MT> 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<typename WalkType>
void correlation_matrix_uniform_sampling(const unsigned int n, const unsigned int num_points, std::string walkname){

Expand All @@ -106,7 +126,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 @@ -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<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};

for(unsigned int n : dimensions){

old_uniform_sampling<BilliardWalk>(n, num_points);
old_uniform_sampling<BilliardWalk>(n, num_points);

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

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

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

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

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

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

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

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

return 0;
}
}
37 changes: 34 additions & 3 deletions include/matrix_operations/EigenvaluesProblems.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -168,8 +168,9 @@ class EigenvaluesProblems<NT, Eigen::Matrix<NT,Eigen::Dynamic,Eigen::Dynamic>, E
Spectra::DenseCholesky<NT> 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<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 +325,10 @@ 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
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 +441,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
1 change: 0 additions & 1 deletion include/random_walks/uniform_billiard_walk.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,6 @@
#define RANDOM_WALKS_UNIFORM_BILLIARD_WALK_HPP

#include <Eigen/Eigen>

#include "convex_bodies/ball.h"
#include "convex_bodies/ballintersectconvex.h"
#include "convex_bodies/hpolytope.h"
Expand Down

0 comments on commit 9a6676a

Please sign in to comment.