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

Question

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

До сих пор я выяснил, что я должен использовать 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 

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

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

---

Update

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

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

Answer 1

На самом деле здесь нет ошибок, но вы сделали ошибку, поставив 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 называется и для этого СНП-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 Вы можете следить за этим вопросом, если вы заинтересованы.

Answer 2

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