Using the GSL to compute eigenvalues

source: http://www.r-bloggers.com/using-the-gsl-to-compute-eigenvalues/?utm_source=feedburner&utm_medium=email&utm_campaign=Feed%3A+RBloggers+%28R+bloggers%29

Two posts showed how to compute eigenvalues using Armadillo andusing Eigen. As we also looked at using the
GNU GSL, this post will show how to conpute eigenvalues using GSL.

As mentioned in the previous GSL post, we instantiate C language pointers suitable for GSL (here the matrix M). Those must be freed manually, as shown before the return statement.

// [[Rcpp::depends(RcppGSL)]]

#include <RcppGSL.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_eigen.h>

// [[Rcpp::export]]
Rcpp::NumericVector getEigenValues(Rcpp::NumericMatrix sM) {

    RcppGSL::matrix<double> M(sM);  // create gsl data structures from SEXP
    int k = M.ncol();
    Rcpp::NumericVector N(k);   // to store results 

    gsl_vector *eigval = gsl_vector_alloc(k);
    gsl_eigen_symm_workspace *w = gsl_eigen_symm_alloc(k);
    gsl_eigen_symm (M, eigval, w);
    gsl_eigen_symm_free (w);

    for (int j = 0; j < k; j++) {
        N[j] = gsl_vector_get(eigval, j);
    }
    M.free();                          // important: GSL wrappers use C structure
    return N;    // return vector  
}

We can illustrate this easily via a random sample matrix.

set.seed(42)
X <- matrix(rnorm(4*4), 4, 4)
Z <- X %*% t(X)

getEigenValues(Z)

[1] 14.2100  2.4099  1.6856  0.3319

In comparison, R gets the same results (in reverse order) and also returns the eigenvectors.

eigen(Z)

$values
[1] 14.2100  2.4099  1.6856  0.3319

$vectors
         [,1]     [,2]    [,3]     [,4]
[1,]  0.69988 -0.55799  0.4458 -0.00627
[2,] -0.06833 -0.08433  0.0157  0.99397
[3,]  0.44100 -0.15334 -0.8838  0.03127
[4,]  0.55769  0.81118  0.1413  0.10493