Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

add spectra routine for minPosLinearEigenvalue_EigenSymSolver #303

Merged
merged 6 commits into from
Jun 6, 2024
Merged

add spectra routine for minPosLinearEigenvalue_EigenSymSolver #303

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
f981787
add spectra routine for minPosLinearEigenvalue_EigenSymSolver
kasryan May 31, 2024
13ac65d
add new function to check for correlational matrices
kasryan Jun 4, 2024
02fec52
check for correlation matrices after sampling function
kasryan Jun 5, 2024
8501164
update tol in epsilon equality check
kasryan Jun 6, 2024
f2e5079
check for matrices with greater dimensions
kasryan Jun 6, 2024
9b7f114
remove <iostream> header
atrayees Jun 6, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
64 changes: 50 additions & 14 deletions examples/correlation_matrices/sampler.cpp
View file
Open in desktop
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
Open in desktop
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
2 changes: 1 addition & 1 deletion include/random_walks/uniform_billiard_walk.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@
#define RANDOM_WALKS_UNIFORM_BILLIARD_WALK_HPP

#include <Eigen/Eigen>

#include <iostream>
TolisChal marked this conversation as resolved.
Show resolved Hide resolved
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Did you maybe forget to delete it?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes I'm sorry for this- I deleted it just now

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