Какой самый быстрый способ вычислить первые два главных компонента в R?

Question

Я использую princomp в R для выполнения PCA. Моя матрица данных огромна (10K x 10K с каждым значением до 4 десятичных знаков). На процессор Xeon 2,27 ГГц требуется ~ 3,5 часа и ~ 6,5 ГБ физической памяти.

Поскольку мне нужны только первые два компонента, есть ли более быстрый способ сделать это?

Обновление:

В дополнение к скорости, есть ли эффективный способ памяти?

Для расчета первых двух компонентов с использованием svd(,2,) требуется ~ 2 часа и ~ 6,3 ГБ физической памяти.

Answer 1

Иногда вы получаете доступ к так называемым «экономическим» разложениям, которые позволяют ограничить количество собственных значений / собственных векторов. Похоже, что eigen() и prcomp() не предлагают этого, но svd() позволяет вам указать максимальное число для вычисления.

На маленьких матрицах усиления кажутся скромными:

R> set.seed(42); N <- 10; M <- matrix(rnorm(N*N), N, N)
R> library(rbenchmark)
R> benchmark(eigen(M), svd(M,2,0), prcomp(M), princomp(M), order="relative")
          test replications elapsed relative user.self sys.self user.child
2 svd(M, 2, 0)          100   0.021  1.00000      0.02        0          0
3    prcomp(M)          100   0.043  2.04762      0.04        0          0
1     eigen(M)          100   0.050  2.38095      0.05        0          0
4  princomp(M)          100   0.065  3.09524      0.06        0          0
R> 

, но коэффициент три относительно princomp() может стоить вашего при реконструкции princomp() из svd(), поскольку svd() позволяет вам остановиться после двух значений.

Answer 2

Пакет 'svd' предоставляет процедуры для усеченного SVD / собственного разложения по алгоритму Ланцоша.Вы можете использовать его для вычисления только первых двух основных компонентов.

Здесь у меня есть:

> library(svd)
> set.seed(42); N <- 1000; M <- matrix(rnorm(N*N), N, N)
> system.time(svd(M, 2, 0))
   user  system elapsed 
  7.355   0.069   7.501 
> system.time(princomp(M))
   user  system elapsed 
  5.985   0.055   6.085 
> system.time(prcomp(M))
   user  system elapsed 
  9.267   0.060   9.368 
> system.time(trlan.svd(M, neig = 2))
   user  system elapsed 
  0.606   0.004   0.614 
> system.time(trlan.svd(M, neig = 20))
   user  system elapsed 
  1.894   0.009   1.910
> system.time(propack.svd(M, neig = 20))
   user  system elapsed 
  1.072   0.011   1.087 

Answer 3

Я попытался реализовать в пакете pcaMethods алгоритм nipals.По умолчанию он рассчитывает первые 2 основных компонента.Оказывается медленнее, чем другие предложенные методы.

set.seed(42); N <- 10; M <- matrix(rnorm(N*N), N, N)
library(pcaMethods)
library(rbenchmark)
m1 <- pca(M, method="nipals", nPcs=2)
benchmark(pca(M, method="nipals"),
          eigen(M), svd(M,2,0), prcomp(M), princomp(M), order="relative")

                       test replications elapsed relative user.self sys.self
3              svd(M, 2, 0)          100    0.02      1.0      0.02        0
2                  eigen(M)          100    0.03      1.5      0.03        0
4                 prcomp(M)          100    0.03      1.5      0.03        0
5               princomp(M)          100    0.05      2.5      0.05        0
1 pca(M, method = "nipals")          100    0.23     11.5      0.24        0

Answer 4

Способ питания может быть тем, что вы хотите. Если вы кодируете его в R, что совсем не сложно, я думаю, вы можете обнаружить, что он не быстрее, чем SVD-подход, предложенный в другом ответе, который использует скомпилированные подпрограммы LAPACK.

Answer 5

Я удивлен, что никто еще не упомянул пакет irlba:

Это даже немного быстрее, чем svd propack.svd, обеспечивает irlba::prcomp_irlba(X, n=2) stats::prcomp -подобный интерфейс для удобства и не требует настройки параметров в следующем тесте для прямоугольных матриц (2: 1) различного размера. Для матриц размером 6000x3000 это в 50 раз быстрее, чем stats::prcomp. Для матриц меньше 100x50 stats::svd все же быстрее, но.

library(microbenchmark)
library(tidyverse)
#install.packages("svd","corpcor","irlba","rsvd")

exprs <- rlang::exprs(
  svd(M, 2, 2)$v,
  prcomp(M)$rotation[,1:2],
  irlba::prcomp_irlba(M, n=2)$rotation,
  irlba::svdr(M, k=2)$v,
  rsvd::rsvd(M, 2)$v,
  svd::propack.svd(M, neig=2, opts=list(maxiter=100))$v,
  corpcor::fast.svd(M)$v[,1:2]
)

set.seed(42)
tibble(N=c(10,30,100,300,1000,3000)) %>%
  group_by(N) %>%
  do({
    M <- scale(matrix(rnorm(.$N*.$N*2), .$N*2, .$N))
    microbenchmark(!!!exprs,
      times=min(100, ceiling(3000/.$N)))%>%
      as_tibble
  }) %>% 
ggplot(aes(x=N, y=time/1E9,color=expr)) +
  geom_jitter(width=0.05) +
  scale_x_log10("matrix size (2N x N)") +
  scale_y_log10("time [s]") +
  stat_summary(fun.y = median, geom="smooth") +
  scale_color_discrete(labels = partial(str_wrap, width=30))

Рандомизированный SVD, предоставленный rsvd, еще быстрее, но, к сожалению, совсем выключен:

set.seed(42)
N <- 1000
M <- scale(matrix(rnorm(N^2*2), N*2, N))
cor(set_colnames(sapply(exprs, function(x) eval(x)[,1]), sapply(exprs, deparse)))

                                                       svd(M, 2, 2)$v prcomp(M)$rotation[, 1:2] irlba::prcomp_irlba(M, n = 2)$rotation irlba::svdr(M, k = 2)$v rsvd::rsvd(M, 2)$v svd::propack.svd(M, neig = 2, opts = list(maxiter = 100))$v corpcor::fast.svd(M)$v[, 1:2]
svd(M, 2, 2)$v                                                   1.0000000                 1.0000000                             -1.0000000               0.9998748           0.286184                                                   1.0000000                     1.0000000
prcomp(M)$rotation[, 1:2]                                        1.0000000                 1.0000000                             -1.0000000               0.9998748           0.286184                                                   1.0000000                     1.0000000
irlba::prcomp_irlba(M, n = 2)$rotation                          -1.0000000                -1.0000000                              1.0000000              -0.9998748          -0.286184                                                  -1.0000000                    -1.0000000
irlba::svdr(M, k = 2)$v                                          0.9998748                 0.9998748                             -0.9998748               1.0000000           0.290397                                                   0.9998748                     0.9998748
rsvd::rsvd(M, 2)$v                                               0.2861840                 0.2861840                             -0.2861840               0.2903970           1.000000                                                   0.2861840                     0.2861840
svd::propack.svd(M, neig = 2, opts = list(maxiter = 100))$v      1.0000000                 1.0000000                             -1.0000000               0.9998748           0.286184                                                   1.0000000                     1.0000000
corpcor::fast.svd(M)$v[, 1:2]                                    1.0000000                 1.0000000                             -1.0000000               0.9998748           0.286184                                                   1.0000000                     1.0000000

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

Answer 6

Пакеты "gmodels" и "corpcor" R поставляются с более быстрой реализацией SVD и PCA. Они работают аналогично основным версиям для маленьких матриц:

> set.seed(42); N <- 10; M <- matrix(rnorm(N*N), N*N, N)
> library("rbenchmark")
> library("gmodels")    
> benchmark(svd(M,2,0), svd(M), gmodels::fast.svd(M), corpcor::fast.svd(M), prcomp(M), gmodels::fast.prcomp(M), princomp(M), order="relative")
                     test replications elapsed relative user.self sys.self user.child sys.child
1            svd(M, 2, 0)          100   0.005      1.0     0.005    0.000          0         0
2                  svd(M)          100   0.006      1.2     0.005    0.000          0         0
3    gmodels::fast.svd(M)          100   0.007      1.4     0.006    0.000          0         0
4    corpcor::fast.svd(M)          100   0.007      1.4     0.007    0.000          0         0
6 gmodels::fast.prcomp(M)          100   0.014      2.8     0.014    0.000          0         0
5               prcomp(M)          100   0.015      3.0     0.014    0.001          0         0
7             princomp(M)          100   0.030      6.0     0.029    0.001          0         0
> 

Однако они обеспечивают более быстрый результат для больших матриц (особенно с множеством строк).

> set.seed(42); N <- 10; M <- matrix(rnorm(N*N), N*N*N, N)
> library("rbenchmark")
> library("gmodels")
> benchmark(svd(M,2,0), svd(M), gmodels::fast.svd(M), corpcor::fast.svd(M), prcomp(M), gmodels::fast.prcomp(M), order="relative")

                     test replications elapsed relative user.self sys.self user.child sys.child
4    corpcor::fast.svd(M)          100   0.029    1.000     0.028    0.001          0         0
3    gmodels::fast.svd(M)          100   0.035    1.207     0.033    0.001          0         0
2                  svd(M)          100   0.037    1.276     0.035    0.002          0         0
1            svd(M, 2, 0)          100   0.039    1.345     0.037    0.001          0         0
5               prcomp(M)          100   0.068    2.345     0.061    0.006          0         0
6 gmodels::fast.prcomp(M)          100   0.068    2.345     0.060    0.007          0         0

Answer 7

вы можете использовать подход нейронной сети, чтобы найти основной компонент.Основное описание дается здесь .. http://www.heikohoffmann.de/htmlthesis/node26.html

Первый главный компонент, y = w1 * x1 + w2 * x2, а второй ортогональный компонент можно вычислить как q = w2 * x1-w1 * x2.

Answer 8

Вы можете написать функцию самостоятельно и остановиться на 2 компонентах.Это не так уж сложно.У меня это где-то валяется, если я найду, я опубликую.