Using the GSL to compute eigenvalues

January 18, 2013

(This article was first published on Rcpp Gallery, and kindly contributed to R-bloggers)

Two posts showed how to compute eigenvalues using Armadillo and using 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();

    RcppGSL::vector<double> 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);

    // we need to invoke wrap() here to help the compiler, but at least we save a loop
    // we also need to create a new Rcpp object as the RcppGSL one needs to be free'd.
    Rcpp::NumericVector res(Rcpp::wrap(ev));;               		// important: GSL wrappers use C structure;

    return res;				// return results vector  

We can illustrate this easily via a random sample matrix.

X <- matrix(rnorm(4*4), 4, 4)
Z <- X %*% t(X)

[1] 14.2100  2.4099  1.6856  0.3319

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

[1] 14.2100  2.4099  1.6856  0.3319

         [,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

To leave a comment for the author, please follow the link and comment on their blog: Rcpp Gallery. offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...

If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.


Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)