Using the GSL to compute eigenvalues

Dirk Eddelbuettel — written Jan 18, 2013 — updated Nov 23, 2015 — source

Two previous posts showed how to compute eigenvalues using Armadillo and using Eigen. As we have 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). Prior to release 0.3.0 of RcppGSL, these had to be freed manually. However, since release 0.3.0 since is now taken care of via the standard C++ mechanism of destructors.

// Tell Rcpp to rely on the RcppGSL package to find GSL library and headers
// [[Rcpp::depends(RcppGSL)]]

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

// [[Rcpp::export]]
RcppGSL::Vector getEigenValues(RcppGSL::Matrix & M) {
    int k = M.ncol();

    RcppGSL::Vector ev(k);  	// instead of gsl_vector_alloc(k);
    gsl_eigen_symm_workspace *w = gsl_eigen_symm_alloc(k);
    gsl_eigen_symm(M, ev, w);
    gsl_eigen_symm_free (w);

    return ev;				// return results 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.210011  2.409920  1.685588  0.331887

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

eigen(Z)

$values
[1] 14.210011  2.409920  1.685588  0.331887

$vectors
           [,1]       [,2]       [,3]        [,4]
[1,]  0.6998802 -0.5579950  0.4458363 -0.00626959
[2,] -0.0683337 -0.0843295  0.0157046  0.99396801
[3,]  0.4409966 -0.1533407 -0.8837594  0.03127147
[4,]  0.5576919  0.8111823  0.1412537  0.10493036

tags: modeling gsl

Extending R with C++ and Fortran — Dirk Eddelbuettel and JBrandon Duck-Mayr
Benchmarking Rcpp code with RcppClock — Zach DeBruine
Simulation Smoother using RcppArmadillo — Tomasz Woźniak

Using the GSL to compute eigenvalues

Related Articles