Skip to content

Commit

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

This reverts commit 9a6676a.
  • Loading branch information
@vissarion
vissarion authored Jun 7, 2024
1 parent 9a6676a commit 1c11035
Show file tree
Hide file tree
Showing 3 changed files with 18 additions and 84 deletions.
64 changes: 14 additions & 50 deletions examples/correlation_matrices/sampler.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -89,26 +89,6 @@ 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 @@ -126,7 +106,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" + std::to_string(n) + ".txt", randPoints);
write_to_file<Point>(walkname + "_matrices.txt", randPoints);
}

template<typename WalkType>
Expand All @@ -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<std::chrono::milliseconds>(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<PointMT>(walkname + "_matrices_MT" + std::to_string(n) + ".txt", randPoints);

write_to_file<PointMT>(walkname + "_matrices_MT.txt", randPoints);
}

int main(int argc, char const *argv[]){
Expand All @@ -177,30 +146,25 @@ int main(int argc, char const *argv[]){
printf("\n");
#endif

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){
unsigned int n = 3, num_points = 5000;

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: 3 additions & 34 deletions include/matrix_operations/EigenvaluesProblems.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -168,9 +168,8 @@ 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, ncv);
geigs(&op, &Bop, 1, 15 < matrixDim ? 15 : matrixDim);

// Initialize and compute
geigs.init();
Expand Down Expand Up @@ -325,10 +324,9 @@ 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 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<NT, Spectra::LARGEST_ALGE, Spectra::DenseSymMatProd<NT>, Spectra::DenseCholesky<NT>, Spectra::GEIGS_CHOLESKY>
geigs(&op, &Bop, 1, ncv);
geigs(&op, &Bop, 1, 15 < matrixDim ? 15 : matrixDim);

// Initialize and compute
geigs.init();
Expand Down Expand Up @@ -441,39 +439,10 @@ 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: 1 addition & 0 deletions include/random_walks/uniform_billiard_walk.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@
#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 1c11035

Please sign in to comment.