Site icon R-bloggers

Armadillo eigenvalues

[This article was first published on Rcpp Gallery, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Today a (slightly confused) question on StackOverflow wondered how to access R’s facilities for eigenvalues calculations from C code.

For this, we need to step back and consider how this is done. In fact, R farms the calculation out to the BLAS. On could possibly access R’s functions—but would then have to wrestle with the data input/output issues which make Rcpp shine in comparison. Also, Rcpp gets us access to Armadillo (via the RcppArmadillo) package and Armadillo’s main focus are exactly the linear algebra calculations and decompositions.

And with facilities that were added to Rcpp in the 0.10.* release series, this effectively becomes a one-liner of code! (Nitpickers will note that there are also one include statement, two attributes declarations and the function name itself.)

#include <RcppArmadillo.h>

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

// [[Rcpp::export]]
arma::vec getEigenValues(arma::mat M) {
    return arma::eig_sym(M);
}

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]
[1,]  0.3319
[2,]  1.6856
[3,]  2.4099
[4,] 14.2100

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

Armadillo has other eigenvector computations too, see its documentation.

To leave a comment for the author, please follow the link and comment on their blog: Rcpp Gallery.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.