An issue regarding symmetry
At this point I had my Rcpp code running, and on testing the function on R, this is what appeared-
result[[1]][[1]]
[,1] [,2] [,3] [,4]
[1,] 1.00000000 0.01231237 -0.16375568 -0.07145539
[2,] 0.01231237 1.0000000 -0.7565504 0.4482418
[3,] -0.16375568 -0.7565504 1.0000000 0.2120664
[4,] -0.07145539 0.4482418 0.2120664 1
(Here is the example code for running the function-)
# Example for sampling correlation matrices
#Import required library
library(volesti)
#Input parameters:
# n: The dimension of the correlation matrix
# num_matrices: The number of matrices to generate (default is 1000)
# walk_length: The length of the random walk (default is 1)
# nburns: The number of burn-in steps for the random walk (default is 0)
# validate: Whether to validate the sampled matrices (default is FALSE)
result <- uniform_sample_correlation_matrices(n = 4, num_matrices = 5000, walk_length = 2, nburns = 0, validate = TRUE)
# result now contains the list of sampled correlation matrices
This is a very undesired output; the problem arises because only the lower triangular part of the matrix is used in C++ computations- and while this improves memory efficiency in the C++ code, it is also setting the upper triangular part to be zero, because it assumes the matrix to be symmetrical.
Correcting the Code
Previously, the function in C++ responsible for the uniform sampling of correlation matrices was carrying out all computations using only the lower triangular part of the matrix. While that is okay, we also needed to fix the upper triangular matrix using the line of code:
` MT final_cor_mat = p.mat + p.mat.transpose() - MT::Identity(n, n); `
That being so, here was my new code for the function:
template
<
typename WalkTypePolicy,
typename PointType,
typename RNGType,
typename MT
>
void uniform_correlation_sampling_MT( const unsigned int &n,
std::list<MT> &randCorMatrices,
const unsigned int &walkL,
const unsigned int &num_points,
unsigned int const& nburns){
using PointList = std::list<PointType>;
PointList randPoints;
CorrelationSpectrahedron_MT<PointType> P(n);
const unsigned int d = P.dimension();
PointType startingPoint(n);
RNGType rng(d);
uniform_sampling<WalkTypePolicy>(randPoints, P, rng, walkL, num_points, startingPoint, nburns);
for(const auto&p : randPoints){
MT final_cor_mat = p.mat + p.mat.transpose() - MT::Identity(n, n);
randCorMatrices.push_back(final_cor_mat);
}
}
One significant change was that previously, the function was taking PointList &randPoints as an input parameter to store the CorreMatrix
Consequently, I applied this change to two other functions as well- ` gaussian_correlation_sampling_MT() ` and ` exponential_correlation_sampling_MT() `.