4

mR で大規模なスパース行列の最初の固有ベクトルを計算しようとしています。ここでeigen()は N > 10 6を意味するため、使用は現実的ではありません。

igraphこれまでのところ、スパース行列を処理できるパッケージの ARPACK を使用する必要があることがわかりました。ただし、非常に単純な (3x3) マトリックスで動作させることはできません。

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 を使用する関数のバグによるものと思われます。

バグを解決するための提案 (パッケージ コードを確認しようとしましたが、複雑すぎて理解できません)、または別の方法で固有ベクトルを計算するための提案を歓迎します。

4

2 に答える 2

4

ここには実際にはバグはありませんが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

詳細に興味がある場合は、対称の代わりに、一般的な非対称固有値ソルバーが呼び出され、そのために 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()非常に遅くなる可能性があります。問題は、反復ごとに C コードから R にコールバックする必要があることです。このスレッドを参照して ください: http://lists.gnu.org/archive/html/igraph-help/2012-02/msg00029.htmlarpack()多くの反復が必要であり、後者は行列の固有構造に関連しています。

R コールバックの代わりに Rcpp を使用してオプションで C コールバックを使用できるかどうかを確認するために、igraph 課題トラッカーで課題を作成しました: https://github.com/igraph/igraph/issues/491 フォローできます興味があればこの問題。

于 2013-06-08T14:50:37.693 に答える