In the fourth and the fifth weeks, I worked on Rvolesti: this is an R interface for the C++ library volesti.
Background
Up till now, my work focused on solving generalized Eigenvalue problems in each step of the random walk- this is required in sampling correlation matrices.
My job was to expose the function that will uniformly sample matrices to the R interface, for the user to make use of this functionality. This functionality enables users in R to generate matrices with uniformly distributed random elements, which can be particularly useful in statistical modeling, and various other applications.
The code was in Rcpp; Rcpp is a powerful R package that allows you to integrate code in R and C++ consistently.
My Code
// [[Rcpp::export]]
Rcpp::List uniform_sample_correlation_matrices(const unsigned int n, const unsigned int num_matrices = 1000,
const unsigned int walk_length=1, const unsigned int nburns = 0, const bool validate = false) {
typedef double NT;
typedef Eigen::Matrix<NT, Eigen::Dynamic, Eigen::Dynamic> MT;
typedef BoostRandomNumberGenerator<boost::mt19937, NT> RNGType;
typedef CorreMatrix<NT> PointMT;
std::list<MT> randCorMatrices;
uniform_correlation_sampling_MT<AcceleratedBilliardWalk, PointMT, RNGType>(n, randCorMatrices, walk_length, num_matrices, nburns);
EigenvaluesProblems<NT, MT, Eigen::Matrix<NT, Eigen::Dynamic, 1>> solver;
if (validate) {
const NT tol = 1e-8;
for (const auto& mat : randCorMatrices) {
if (!solver.is_correlation_matrix(mat, tol)) {
throw Rcpp::exception("Invalid correlation matrix");
}
}
}
Rcpp::List rcpp_sampled_matrices(randCorMatrices.size());
int index = 0;
for (const auto& mat : randCorMatrices) {
rcpp_sampled_matrices[index++] = Rcpp::wrap(mat);
}
Rcpp::List result = Rcpp::List::create(
Rcpp::Named("sampled_matrices") = rcpp_sampled_matrices
);
return result;
}
This Rcpp function uniformly samples correlation matrices. It has the following input parameters:
- n; the dimension of the correlation matrix
- num_matrices; the number of matrices to generate
- walk_length; the length of the random walk
- nburns; the number of burn-in steps for the random walk
- validate; this variable will indicate whether or not to validate the sample matrices. If true, then the function checks if all sampled matrices are indeed correlation matrices.
The function returns a list of sampled matrices.
The C++ function, uniform_correlation_sampling_MT() is called to generate the random correlation matrices, which are stored in the list randCorMatrices. This is followed by an optional validation loop, which will throw an error if any matrix in randCorMatrices is not a correlation matrix.
Now, it is important that the data structures used in the C++ code are all compatible with the R environment for proper integration of the Rcpp code. This is why the Eigen matrices are converted to Rcpp::List here- this ensures that the sampled correlation matrices generated by the C++ function are returned in a format that R users can readily use and manipulate. Rcpp::wrap() handles the conversion.
The ` // [[Rcpp::export]] ` attribute specifies that this function is to be exported to R. So, on compilation of this function-
` Rcpp::compileAttributes() ` ,
the RcppExports.cpp and RcppExports.R files are generated. Now, the code becomes usable in the R interface when the volesti package is built again in R.
The link to my pull request is here.
