Вычисление собственных значений в кодах C-внутри-R

Если вы готовы использовать Rcpp и связанные с ним пакеты, существующие примеры для fastLm() покажут вам, как это сделать с помощью

• Eigen через RcppEigen
• Армадилло через Rcpp
• GSL через RcppGSL

где последние два будут использовать тот же BLAS, что и R, а в Eigen есть что-то внутреннее, что часто может быть быстрее. Все пакеты реализуют lm() с использованием декомпозиции, предлагаемой языком (часто с использованием решения или связанного, но переключиться на SVD несложно, если у вас есть инструментальная цепочка, работающая на вас).

Изменить: вот явный пример. Используйте следующий код C ++:

#include <RcppArmadillo.h>   

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

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

сохраните в файл «eigenEx.cpp». Тогда все, что нужно, это этот код R:

library(Rcpp)               ## recent version for sourceCpp()
sourceCpp("eigenEx.cpp")    ## converts source file into getEigen() we can call

чтобы мы могли запустить это:

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

getEigen(Z)                 ## calls function created above

и я получаю тот же собственный вектор, что и у R. Это действительно не намного проще.

Это также позволяет Armadillo выбирать, какой метод использовать для разложения Eigen, и, как намекнул Дэвид, это что-то более быстрое, чем полномасштабная SVD (см. Документацию Armadillo).

Вы можете использовать одну из многих доступных библиотек линейной алгебры (Lapack). Здесь есть ссылка, объясняющая, как получить библиотеки Lapack для Windows. GOTOBLAS и ACML бесплатны. MKL также бесплатен для некоммерческого использования. После того, как вы установили библиотеки, вам нужна функция sgesvd (для float) или dgesvd (для double).

Вот пара примеров от Intel о том, как использовать gesvd.

1. Row Major
2. Col ​​Major

Если вы используете Linux, ознакомьтесь с GNU SL и Eigen. Эти библиотеки обычно можно установить из диспетчера пакетов дистрибутива.