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

Feature/spectra correlations #306

Merged
merged 13 commits into from
Jun 17, 2024
Merged

Feature/spectra correlations #306

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
12359ad
update the ncv values
atrayees Jun 10, 2024
742df12
include spectra routine
atrayees Jun 10, 2024
6e0a6bc
function to check for a correlation matrix
atrayees Jun 10, 2024
b36ea3d
run sampling for matrices with greater dimensions
atrayees Jun 10, 2024
93e88c8
remove blank space lines
atrayees Jun 10, 2024
e7e8252
fix indentations
atrayees Jun 12, 2024
322f794
fixed indentation and code style
atrayees Jun 12, 2024
56ea104
updated tol as a function parameter
atrayees Jun 12, 2024
594e1d0
fix indentations
atrayees Jun 12, 2024
91b10d4
update comments about the ncv parameter
atrayees Jun 12, 2024
f872b12
Update comment regarding computing of the minimum positive eigenvalue…
atrayees Jun 12, 2024
2b52486
reuse function isPositiveSemidefinite()
atrayees Jun 14, 2024
8128efa
fix indentation
atrayees Jun 14, 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
59 changes: 46 additions & 13 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){
vissarion marked this conversation as resolved.
Show resolved Hide resolved

const double tol = 1e-8;
vissarion marked this conversation as resolved.
Show resolved Hide resolved

//check if all the diagonal elements are ones
for(int i=0 ; i<matrix.rows() ; i++){
vissarion marked this conversation as resolved.
Show resolved Hide resolved
if(std::abs(matrix(i, i)-1.0) > tol){
return false;
}
}

//check if the matrix is positive definite
Copy link
Contributor

Choose a reason for hiding this comment

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

Why not using the function isPositiveSemidefinite. Note: that function is not using tolerance.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I have updated the function and made use of the isPositiveSemidefinite function

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;
vissarion marked this conversation as resolved.
Show resolved Hide resolved
}

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 @@ -126,7 +146,15 @@ void correlation_matrix_uniform_sampling_MT(const unsigned int n, const unsigned
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++;
Copy link
Member

Choose a reason for hiding this comment

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

Please fix indentation here

Copy link
Contributor Author

Choose a reason for hiding this comment

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

done

}
}
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 +174,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};

old_uniform_sampling<BilliardWalk>(n, num_points);
for(unsigned int n : dimensions){

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

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

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

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

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

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

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

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

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
Copy link
Member

Choose a reason for hiding this comment

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

"... computing the minimum positive eigenvalue by computing the largest eigenvalue of the inverse Generalized Eigenvalue Problem."

Copy link
Contributor Author

Choose a reason for hiding this comment

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

i updated the comment.

int ncv = std::min(std::max(10, matrixDim/20), matrixDim);
Copy link
Contributor

Choose a reason for hiding this comment

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

Could you please comment on where those numbers come from? 10, 20, 15 ?

Copy link
Member

Choose a reason for hiding this comment

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

This is an empirical value I gave to @atrayees. We left as a task after the R interface to tune this implementation (not only ncv value).

Copy link
Member

Choose a reason for hiding this comment

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

Thanks @TolisChal for the explanations. @atrayees could you please add a comment with this explanation and mention that as a TODO?

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, done

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();
vissarion marked this conversation as resolved.
Show resolved Hide resolved
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);
vissarion marked this conversation as resolved.
Show resolved Hide resolved
}

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
Loading