跳到内容

ARPACK 特征向量计算

源码：R/centrality.R
arpack.Rd

用于计算稀疏矩阵特征向量的 ARPACK 库接口

用法

arpack_defaults()

arpack(
  func,
  extra = NULL,
  sym = FALSE,
  options = arpack_defaults(),
  env = parent.frame(),
  complex = !sym
)

参数

func

执行矩阵-向量乘法的函数。 ARPACK 需要用户执行这些操作。该函数将向量 \(x\) 作为第一个参数，并且应该返回 \(Ax\)，其中 \(A\) 是“输入矩阵”。（输入矩阵永远不会被显式给出。）第二个参数是 extra。

extra

要提供给 func 的额外参数。

sym

逻辑标量，输入矩阵是否对称。如果矩阵对称，请始终在此处提供 TRUE，因为它可以加快计算速度。

options

ARPACK 的选项，一个命名列表，用于覆盖一些默认选项值。请参见下面的详细信息。

env

将评估 func 的环境。

complex

是否将 ARPACK 返回的特征向量转换为 R 复向量。默认情况下，对于对称问题（这些问题只有实特征向量/值）不会这样做，而只会对非对称问题这样做。如果您有一个非对称问题，但您确定结果将是实数，那么请在此处提供 FALSE。

值

一个包含以下成员的命名列表

values

数值向量，所需的特征值。

vectors

数值矩阵，所需的特征向量作为列。如果 complex=TRUE（非对称问题的默认设置），则矩阵是复数。

options

一个命名列表，其中包含提供的 options 和一些关于执行的计算的信息，包括 ARPACK 退出代码。请参阅上面的详细信息。

详细信息

ARPACK 是一个用于解决大规模特征值问题的库。该软件包旨在计算一般 \(n\) 乘以 \(n\) 矩阵 \(A\) 的几个特征值和相应的特征向量。它最适合于大型稀疏或结构化矩阵 \(A\)，其中结构化意味着矩阵-向量乘积 w <- Av 需要阶 \(n\) 而不是通常的阶 \(n^2\) 浮点运算。

此函数是 ARPACK 的接口。 igraph 不包含所有 ARPACK 例程，仅包含使用双精度实数处理对称和非对称特征值问题的例程。

ARPACK 中的特征值计算（在最简单的情况下）涉及计算 \(Av\) 乘积，其中 \(A\) 是我们使用的矩阵，\(v\) 是任意向量。在 fun 参数中提供的函数应执行此乘积。如果可以有效地完成乘积，例如，如果矩阵是稀疏的，那么 arpack() 通常能够非常快速地计算特征值。

options 参数指定要执行的计算类型。它是一个包含以下成员的列表，它们直接对应于 ARPACK 参数。在输入时，它具有以下字段

bmat

字符常量，可能的值：‘I’，标准特征值问题，\(Ax=\lambda x\)；和 ‘G’，广义特征值问题，\(Ax=\lambda B x\)。目前仅支持 ‘I’。

n

数值标量。特征值问题的维度。仅当您直接调用 arpack() 时才需要设置此参数。（即，对于 eigen_centrality(), page_rank() 等，不需要此参数。）

which

指定要计算的特征值/向量，包含正好两个字符的字符常量。对称输入矩阵的可能值

"LA"

计算 nev 个最大（代数）特征值。

"SA"

计算 nev 个最小（代数）特征值。

"LM"

计算 nev 个最大（幅度）特征值。

"SM"

计算 nev 个最小（幅度）特征值。

"BE"

计算 nev 个特征值，每个端点一半。当 nev 为奇数时，从高端计算的比从低端计算的多一个。

非对称输入矩阵的可能值

"LM"

计算 nev 个最大幅度的特征值。

"SM"

计算 nev 个最小幅度的特征值。

"LR"

计算 nev 个最大实部的特征值。

"SR"

计算 nev 个最小实部的特征值。

"LI"

计算 nev 个最大虚部的特征值。

"SI"

计算 nev 个最小虚部的特征值。

此参数有时会被各种函数覆盖，例如，page_rank() 始终设置为 ‘LM’。

nev

数值标量。要计算的特征值的数量。

tol

数值标量。停止准则：如果 Ritz 值的误差小于其估计值的 tol 倍，则认为其相对精度是可以接受的。如果将其设置为零，则使用机器精度。

ncv

要生成的 Lanczos 向量的数量。

ldv

数值标量。在当前实现中，应将其设置为零。

ishift

零或一。如果为零，则移位由用户通过反向通信提供。如果为一，则相对于缩减的三对角矩阵 \(T\) 的精确移位。请始终将其设置为一。

maxiter

允许的最大 Arnoldi 更新迭代次数。

nb

要在递归中使用的块大小。请始终将此值保留为默认值，即一。

mode

要解决的特征值问题的类型。如果输入矩阵是对称的，则可能的值

1

\(Ax=\lambda x\)，\(A\) 是对称的。

2

\(Ax=\lambda Mx\)，\(A\) 是对称的，\(M\) 是对称正定的。

3

\(Kx=\lambda Mx\)，\(K\) 是对称的，\(M\) 是对称正半定的。

4

\(Kx=\lambda KGx\)，\(K\) 是对称正半定的，\(KG\) 是对称不定的。

5

\(Ax=\lambda Mx\)，\(A\) 是对称的，\(M\) 是对称正半定的。（Cayley 变换模式。）

请注意，仅测试了 mode==1，其他值可能无法正常工作。如果输入矩阵不对称，则可能的值

1

\(Ax=\lambda x\)。

2

\(Ax=\lambda Mx\)，\(M\) 是对称正定的。

3

\(Ax=\lambda Mx\)，\(M\) 是对称半正定的。

4

\(Ax=\lambda Mx\)，\(M\) 是对称半正定的。

请注意，仅测试了 mode==1，其他值可能无法正常工作。

start

当前未使用。稍后将用于设置起始向量。

sigma

当前未使用。

sigmai

当前未使用。

：在输出时，添加以下附加字段

info

ARPACK 的错误标志。可能的值

0

正常退出。

1

达到最大迭代次数。

3

在隐式重启的 Arnoldi 迭代循环中无法应用移位。一种可能性是相对于 nev 增加 ncv 的大小。

ARPACK 可以返回比这些更多的错误条件，但它们会被转换为常规的 igraph 错误。

iter

进行的 Arnoldi 迭代次数。

nconv

“收敛”的 Ritz 值的数量。这表示满足收敛标准的 Ritz 值的数量。

numop

矩阵-向量乘法的总数。

numopb

当前未使用。

numreo

重新正交化的总步数。

有关更多详细信息，请参见 ARPACK 文档。

参考文献

D.C. Sorensen, Implicit Application of Polynomial Filters in a k-Step Arnoldi Method. SIAM J. Matr. Anal. Apps., 13 (1992), pp 357-385.

R.B. Lehoucq, Analysis and Implementation of an Implicitly Restarted Arnoldi Iteration. Rice University Technical Report TR95-13, Department of Computational and Applied Mathematics.

B.N. Parlett & Y. Saad, Complex Shift and Invert Strategies for Real Matrices. Linear Algebra and its Applications, vol 88/89, pp 575-595, (1987).

参见

eigen_centrality()、page_rank()、hub_score()、cluster_leading_eigen() 是 igraph 中一些使用 ARPACK 的函数。

作者

Rich Lehoucq、Kristi Maschhoff、Danny Sorensen、Chao Yang 为 ARPACK，Gabor Csardi csardi.gabor@gmail.com 为 R 接口。

示例


# Identity matrix
f <- function(x, extra = NULL) x
arpack(f, options = list(n = 10, nev = 2, ncv = 4), sym = TRUE)
#> $values
#> [1] 1 1
#> 
#> $vectors
#>              [,1]        [,2]
#>  [1,]  0.47899559 -0.22165241
#>  [2,]  0.02563875 -0.05909943
#>  [3,]  0.47539054 -0.19000747
#>  [4,] -0.52783666 -0.67727036
#>  [5,] -0.27844498 -0.21112117
#>  [6,]  0.14568395  0.04884236
#>  [7,]  0.07307229  0.02239947
#>  [8,] -0.37430168  0.63030259
#>  [9,] -0.02487496  0.05239496
#> [10,]  0.14311774  0.07136832
#> 
#> $options
#> $options$bmat
#> [1] "I"
#> 
#> $options$n
#> [1] 10
#> 
#> $options$which
#> [1] "XX"
#> 
#> $options$nev
#> [1] 2
#> 
#> $options$tol
#> [1] 0
#> 
#> $options$ncv
#> [1] 4
#> 
#> $options$ldv
#> [1] 0
#> 
#> $options$ishift
#> [1] 1
#> 
#> $options$maxiter
#> [1] 3000
#> 
#> $options$nb
#> [1] 1
#> 
#> $options$mode
#> [1] 1
#> 
#> $options$start
#> [1] 0
#> 
#> $options$sigma
#> [1] 0
#> 
#> $options$sigmai
#> [1] 0
#> 
#> $options$info
#> [1] 0
#> 
#> $options$iter
#> [1] 1
#> 
#> $options$nconv
#> [1] 2
#> 
#> $options$numop
#> [1] 4
#> 
#> $options$numopb
#> [1] 0
#> 
#> $options$numreo
#> [1] 4
#> 
#> 

# Graph laplacian of a star graph (undirected), n>=2
# Note that this is a linear operation
f <- function(x, extra = NULL) {
  y <- x
  y[1] <- (length(x) - 1) * x[1] - sum(x[-1])
  for (i in 2:length(x)) {
    y[i] <- x[i] - x[1]
  }
  y
}

arpack(f, options = list(n = 10, nev = 1, ncv = 3), sym = TRUE)
#> $values
#> [1] 10
#> 
#> $vectors
#>  [1]  0.9486833 -0.1054093 -0.1054093 -0.1054093 -0.1054093 -0.1054093
#>  [7] -0.1054093 -0.1054093 -0.1054093 -0.1054093
#> 
#> $options
#> $options$bmat
#> [1] "I"
#> 
#> $options$n
#> [1] 10
#> 
#> $options$which
#> [1] "XX"
#> 
#> $options$nev
#> [1] 1
#> 
#> $options$tol
#> [1] 0
#> 
#> $options$ncv
#> [1] 3
#> 
#> $options$ldv
#> [1] 0
#> 
#> $options$ishift
#> [1] 1
#> 
#> $options$maxiter
#> [1] 3000
#> 
#> $options$nb
#> [1] 1
#> 
#> $options$mode
#> [1] 1
#> 
#> $options$start
#> [1] 0
#> 
#> $options$sigma
#> [1] 0
#> 
#> $options$sigmai
#> [1] 0
#> 
#> $options$info
#> [1] 0
#> 
#> $options$iter
#> [1] 1
#> 
#> $options$nconv
#> [1] 1
#> 
#> $options$numop
#> [1] 3
#> 
#> $options$numopb
#> [1] 0
#> 
#> $options$numreo
#> [1] 2
#> 
#> 

# double check
eigen(laplacian_matrix(make_star(10, mode = "undirected")))
#> eigen() decomposition
#> $values
#>  [1] 1.000000e+01 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#>  [6] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 3.552714e-15
#> 
#> $vectors
#>             [,1]          [,2]          [,3]        [,4]          [,5]
#>  [1,]  0.9486833  0.000000e+00  0.000000e+00  0.00000000  0.000000e+00
#>  [2,] -0.1054093 -3.700743e-17 -3.700743e-17  0.00000000 -3.700743e-17
#>  [3,] -0.1054093 -3.768958e-02  1.005977e-01  0.02684135  4.687454e-02
#>  [4,] -0.1054093  3.537753e-01 -1.012990e-01 -0.08701190 -6.664397e-02
#>  [5,] -0.1054093 -4.526060e-01 -4.556778e-01 -0.51370736 -3.394283e-01
#>  [6,] -0.1054093 -3.384853e-01 -1.777216e-03  0.02479089  8.481572e-01
#>  [7,] -0.1054093 -8.057222e-02 -1.154307e-02  0.01661573 -1.545983e-01
#>  [8,] -0.1054093  5.759880e-01 -4.900096e-01  0.30166795  6.743751e-02
#>  [9,] -0.1054093  3.180750e-01  6.685471e-01 -0.43654213 -4.340406e-02
#> [10,] -0.1054093 -3.384853e-01  2.911619e-01  0.66734547 -3.583947e-01
#>                [,6]          [,7]          [,8]       [,9]      [,10]
#>  [1,]  0.000000e+00  0.000000e+00  0.000000e+00  0.0000000 -0.3162278
#>  [2,]  3.238150e-17  7.401487e-17  3.700743e-17  0.9428090 -0.3162278
#>  [3,] -1.979763e-02  9.262957e-01  4.611153e-02 -0.1178511 -0.3162278
#>  [4,]  8.482088e-01 -8.673487e-02  2.424907e-02 -0.1178511 -0.3162278
#>  [5,] -1.084886e-01 -8.673487e-02  2.531916e-01 -0.1178511 -0.3162278
#>  [6,]  4.196770e-02 -1.931187e-01  3.719095e-02 -0.1178511 -0.3162278
#>  [7,] -1.084886e-01 -8.673487e-02 -9.082431e-01 -0.1178511 -0.3162278
#>  [8,] -4.229128e-01 -8.673487e-02  1.455981e-01 -0.1178511 -0.3162278
#>  [9,] -2.724565e-01 -1.931187e-01  1.513033e-01 -0.1178511 -0.3162278
#> [10,]  4.196770e-02 -1.931187e-01  2.505985e-01 -0.1178511 -0.3162278
#> 

## First three eigenvalues of the adjacency matrix of a graph
## We need the 'Matrix' package for this
library("Matrix")
set.seed(42)
g <- sample_gnp(1000, 5 / 1000)
M <- as_adjacency_matrix(g, sparse = TRUE)
f2 <- function(x, extra = NULL) {
  cat(".")
  as.vector(M %*% x)
}
baev <- arpack(
  f2,
  sym = TRUE,
  options = list(
    n = vcount(g),
    nev = 3,
    ncv = 8,
    which = "LM",
    maxiter = 2000
  )
)
#> .....................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................