Eigendecomposition using 'faer' with R (via 'rextendr')

r-faer-eigendecomposition.R

library(rextendr)

rust_source(
  code         = readLines("r-faer-eigendecomposition.rs"),
  dependencies = list("faer" = "0.22.6"),
  features     = "faer",
  profile      = "release"
)

X <- matrix(rnorm(n = 25), 5, 5)

eigen_faer(X)

r-faer-eigendecomposition.rs

use extendr_api::prelude::*;
use ::faer::{Mat, MatRef, Side};
use ::faer::diag::DiagRef;
use ::faer::{Par, set_global_parallelism};
use core::num::NonZeroUsize;

#[extendr]
fn eigen_faer(x: RMatrix<f64>) -> List {

    set_global_parallelism(Par::Rayon(NonZeroUsize::new(12).expect("Value must be nonzero")));
  
    // Convert to a faer matrix
    let A: Mat<f64> = Mat::from_fn(x.nrows(), x.ncols(), |row, col| x[[row, col]]);;

    // Eigendecomposition
    let evd = A.self_adjoint_eigen(Side::Lower).expect("Eigen problem!");

    // Extract matrix of eigenvectors and column matrix of eigenvalues
    let eigenvectors: &MatRef<f64> = &evd.U();
    let eigenvalues: &DiagRef<f64> = &evd.S();
    let dimension: usize = eigenvectors.nrows();
    
    // Convert to types that can be passed to R
    let eigenvectors: RMatrix<f64> = RMatrix::new_matrix(
      dimension, dimension, 
      |row, col| eigenvectors[(row, col)]
    );
    let eigenvalues: Vec<f64> = (0..dimension)
      .map(|i| eigenvalues[i])
      .collect();
    
    // Combine into an R list object
    list!(vectors = eigenvectors, values = eigenvalues)
}