I started coding for the first week of Google Summer of Code from May 27, 2024. My primary focus on this was on GeomScale’s volesti repository on github- mainly on improving a subroutine that computes the smallest possible eigenvalue for a given pair of symmetric matrices.
Background
When working on complex data, you might need to solve a GEP to find the smallest positive eigenvalue- this function does exactly that.
NT minPosLinearEigenvalue_EigenSymSolver(const MT & A, const MT & B, VT &eigvec) const
Say you have a symmetric positive definite matrix A and a symmetric matrix B. You can determine the smallest positive eigenvalue and its corresponding eigenvector efficiently enough using Eigen’s Eigen::GeneralizedSelfAdjointEigenSolver!
Here are the key parts of the code broken down to explain how this works:
- First, you initialize the Eigenvalue Solver.
NT lambdaMinPositive = NT(0); Eigen::GeneralizedSelfAdjointEigenSolver<MT> ges(B, A); - Next, compute and then retrieve the eigenvalues and eigenvectors.
lambdaMinPositive = 1 / ges.eigenvalues().reverse()[0]; eigvec = ges.eigenvectors().reverse().col(0).reverse();
Now where in the project will this come in handy?
volesti implements various algorithms for high-dimensional sampling and volume approximation as well as functions for copula estimation in financial modeling and metabolic network analysis. A part of this lies in sampling uniformly from the set of correlation matrices.
Now, sampling from correlation matrices requires sticking to certain constraints. These constraints are geometrically represented in the form of boundaries of a spectrahedron- the spectrahedron represents the set of all positive semidefinite matrices. By using a ray, you can determine where a point r (representing a matrix) and a direction v intersect the boundary of the spectrahedron. (This ensures that the sampled correlation matrices are always under the constraints.)
Here is where the required application comes in! The minimal positive eigenvalue of a matrix equation gives you the intersection point.
My Code
The pre-existing subroutine used only Eigen routines- my job was to add an implementation of the same function, but using Spectra subroutines, to solve the GEP in each step of the random walk.
This is primarily because while Eigen and Spectra can both be used for computation of Generalized Eigen problems, Spectra is lightweight and is specially optimized for this task- it may provide a better performance for large scale GEP computations.
