Вычисление собственных векторов разреженной матрицы в R

Я пытаюсь вычислить m первые собственные векторы большой разреженной матрицы в R. Использование eigen() не является реалистичным, потому что большое означает N> 106 здесь.

До сих пор я понял, что я должен использовать ARPACK из igraph пакет, который может иметь дело с разреженными матрицами. Однако я не могу заставить его работать на очень простой (3х3) матрице:

library(Matrix)
library(igraph)

TestDiag <- Diagonal(3, 3:1)
TestMatrix <- t(sparseMatrix(i = c(1, 1, 2, 2, 3), j = c(1, 2, 1, 2, 3), x = c(3/5, 4/5, -4/5, 3/5, 1)))
TestMultipliedMatrix <- t(TestMatrix) %*% TestDiag %*% TestMatrix

А затем с помощью кода, приведенного в примере помощи arpack() функция для извлечения 2 первых собственных векторов:

func <- function(x, extra=NULL) { as.vector(TestMultipliedMatrix %*% x) } 
arpack(func, options=list(n = 3, nev = 2, ncv = 3, sym=TRUE, which="LM", maxiter=200), complex = FALSE)

Я получаю сообщение об ошибке:

Error in arpack(func, options = list(n = 3, nev = 2, ncv = 3, sym = TRUE,  :
  At arpack.c:1156 : ARPACK error, NCV must be greater than NEV and less than or equal to N

Я не понимаю эту ошибку, так как ncv (3) здесь больше, чем nev (2), и равно N (3).

Я делаю какую-то глупую ошибку или есть лучший способ для вычисления собственных векторов разреженной матрицы в R?

Обновить

Эта ошибка, по-видимому, связана с ошибкой в arpack() функция с прописными / строчными NCV и NEV.

Любые предложения по устранению ошибки (я попытался взглянуть на код пакета, но он слишком сложен для меня) или вычислить собственные векторы другим способом, приветствуются.

Источник

2 ответа

Решение

На самом деле здесь нет ошибок, но вы ошиблись, поставив sym=TRUE в список опций ARPACK, но sym является аргументом arpack() функция. Т.е. правильный вызов это:

ev <- arpack(func, options=list(n=3, nev=2, ncv=3, which="LM", maxiter=200), 
             sym=TRUE, complex = FALSE)
ev$values
# [1] 3 2
ev$vectors
#               [,1]          [,2]
# [1,] -6.000000e-01 -8.000000e-01
# [2,]  8.000000e-01 -6.000000e-01
# [3,]  2.220446e-16 -9.714451e-17

Если вас интересуют детали, то происходит то, что вместо симметричного вызывается общий несимметричный eigensolver, и для этого также требуется NCV-NEV >= 2. Из источника ARPACK (dnaupd.f):

...
c          NOTE: 2 <= NCV-NEV in order that complex conjugate pairs of Ritz 
c          values are kept together. (See remark 4 below)
...

Еще несколько комментариев, только слабо связанных с вашим вопросом. arpack() может быть довольно медленным Проблема в том, что вам нужно перезванивать в R из кода C на каждой итерации. Смотрите эту ветку: http://lists.gnu.org/archive/html/igraph-help/2012-02/msg00029.html Суть в том, что arpack() Помогает только в том случае, если обратный вызов вашего продукта матрицы-вектора быстрый и вам не нужно много итераций, причем последний связан с собственной структурой матрицы.

Я создал проблему в трекере проблем igraph, чтобы посмотреть, возможно ли дополнительно использовать обратный вызов C, используя Rcpp вместо обратного вызова R: https://github.com/igraph/igraph/issues/491 этот вопрос, если вы заинтересованы.

Источник

Ну, это может быть немного раздражающим, но это работает, когда вы меняете nev=2, ncv=3 в NEV=3, NCV=2, R чувствителен к регистру, что могло вызвать проблему.

Источник
Другие вопросы по тегам