在R中创建Vandermonde矩阵

zaqlnxep  于 2023-07-31  发布在  其他
关注(0)|答案(3)|浏览(130)

我是R的新手,还在学习。我想要m-by-n范德蒙矩阵
x1c 0d1x的数据
我知道可以通过for循环将值赋给矩阵中的相应索引,但当mn很大时,这似乎效率低下。我需要一些建议,有一个更简单,更有效的方法来生成范德蒙矩阵?提前感谢您的帮助!

qyyhg6bp

qyyhg6bp1#

使用outer如下:

n <- 6; alpha <- 1:5 # test data

outer(alpha, seq(0, n-1), `^`)
##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    1    1    1    1    1    1
## [2,]    1    2    4    8   16   32
## [3,]    1    3    9   27   81  243
## [4,]    1    4   16   64  256 1024
## [5,]    1    5   25  125  625 3125

字符串

yptwkmov

yptwkmov2#

有一种基本R方法,它似乎比ThomasIsCoding和G. Grothendieck,这是直接利用^自然矢量化的优势创建矩阵。

vndrMat <- function(x, n) {
  matrix(rep(x, each = n) ^ (seq_len(n) - 1), ncol = n, byrow = TRUE)
}

字符串
下面比较了四个例子:

n <- 6
alpha <- 1:5

vndrMat <- function(x, n) matrix(rep(x, each = n) ^ (seq_len(n) - 1), ncol = n, byrow = TRUE)

out <- function(x, n) outer(x, seq_len(n) - 1, `^`)

vander <- function(alpha,n) t(sapply(alpha, function(k) c(1, cumprod(rep(k, n - 1)))))

library(microbenchmark)
microbenchmark(vndrMat(alpha, n),
               matrixcalc::vandermonde.matrix(alpha, n),
               out(alpha, n),
               vander(alpha, n),
               check = 'identical',
               control = list(order = 'block'),
               times = 1000L)

Unit: microseconds
                                     expr    min     lq      mean median     uq     max neval cld
                        vndrMat(alpha, n)  4.000  4.200  4.343002  4.201  4.302  20.601  1000 a  
 matrixcalc::vandermonde.matrix(alpha, n)  6.301  6.601  6.843913  6.701  6.900  22.101  1000  b 
                            out(alpha, n)  6.000  6.500  6.882389  6.701  6.901  37.702  1000  b 
                         vander(alpha, n) 23.700 24.401 25.566277 24.900 25.501 107.801  1000   c

用@ThomasIsCoding的vander2更新2023-07-10。

托马斯的版本现在是最快的,在统计上与vndrMat并列,但在100,000次迭代中要快一点。

n <- 6
alpha <- 1:5

vndrMat <- function(x, n) matrix(rep(x, each = n) ^ (seq_len(n) - 1), ncol = n, byrow = TRUE)

out <- function(x, n) outer(x, seq_len(n) - 1, `^`)

vander <- function(x, n) t(sapply(x, function(k) c(1, cumprod(rep(k, n - 1)))))

vander2 <- function(x, n) (m <- matrix(x, length(x), n)) ^ (col(m) - 1)

library(microbenchmark)
vTest <- microbenchmark(vndrMat(alpha, n),
                        matrixcalc::vandermonde.matrix(alpha, n),
                        out(alpha, n),
                        vander(alpha, n),
                        vander2(alpha, n),
                        check = 'identical',
                        control = list(order = 'block'),
                        times = 100000L)

print(vTest, order = 'median')

Unit: microseconds
                                     expr  min   lq      mean median   uq     max neval cld
                        vander2(alpha, n)  2.9  3.1  3.550509    3.2  3.3  4599.8 1e+05 a  
                        vndrMat(alpha, n)  2.9  3.2  3.788426    3.3  3.4  6493.3 1e+05 a  
 matrixcalc::vandermonde.matrix(alpha, n)  4.9  5.4  6.183432    5.5  5.8  4544.5 1e+05  b 
                            out(alpha, n)  5.0  5.5  6.343460    5.7  5.9  5766.2 1e+05  b 
                         vander(alpha, n) 22.6 23.8 27.500120   24.4 25.2 95016.6 1e+05   c

> sessionInfo()
R Under development (unstable) (2023-06-22 r84597 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 10 x64 (build 19045)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.utf8  LC_CTYPE=English_United States.utf8   
[3] LC_MONETARY=English_United States.utf8 LC_NUMERIC=C                          
[5] LC_TIME=English_United States.utf8    

time zone: Asia/Jerusalem
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] microbenchmark_1.4.10

loaded via a namespace (and not attached):
 [1] MASS_7.3-60      TH.data_1.1-2    zoo_1.8-12       compiler_4.4.0   Matrix_1.6-0    
 [6] multcomp_1.4-25  sandwich_3.0-2   tools_4.4.0      survival_3.5-5   mvtnorm_1.2-2   
[11] codetools_0.2-19 splines_4.4.0    grid_4.4.0       matrixcalc_1.0-6 lattice_0.21-8

7kqas0il

7kqas0il3#

  • 底座R:sapply + cumprod

一个基本的R解决方案是定义自定义函数vander,其中使用sapply + cumprod

vander1 <- function(alpha,n) t(sapply(alpha, function(k) c(1,cumprod(rep(k,n-1)))))
vm1 <- vander1(alpha,n)

字符串

  • 另一个碱基R:matrix + col
vander2 <- function(alpha, n) (m <- matrix(alpha, length(alpha), n))^(col(m) - 1)
vm2 <- vander2(alpha,n)

  • 从包matrixcalcvandermonde.matrix可以使它
vm3 <- matrixcalc::vandermonde.matrix(alpha,n)

示例

给定alphan,如下所示

alpha <- 1:4
n <- 5


你就会得到

> vm1
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    1    1    1    1
[2,]    1    2    4    8   16
[3,]    1    3    9   27   81
[4,]    1    4   16   64  256
> vm2
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    1    1    1    1
[2,]    1    2    4    8   16
[3,]    1    3    9   27   81
[4,]    1    4   16   64  256
> vm3
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    1    1    1    1
[2,]    1    2    4    8   16
[3,]    1    3    9   27   81
[4,]    1    4   16   64  256

相关问题