我正在尝试对R函数进行编程,以创建一个称为Ui的特殊矩阵(第5页):http : //joshuachan.org/papers/Chan-Jeliazkov-2009.pdf
到目前为止,我已经开发了此功能,我怀疑可以对其进行改进(即提高效率):
create_Uj <- function(uj) {
q <- length(uj)
if (q == 1) return(0)
nr <- q
nc <- q*(q-1)/2
Uj <- matrix(0, nrow = nr, ncol = nc)
for (kk in 2:nr) {
uj_sub <- uj[1:(kk-1)]
Uj[kk, 1:(kk*(kk-1)/2)] <- c(rep(0, (kk-1)*(kk-2)/2), uj[1:(kk-1)])
}
-Uj
}
create_Uj(1:4)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0 0 0 0 0 0
[2,] -1 0 0 0 0 0
[3,] 0 -1 -2 0 0 0
[4,] 0 0 0 -1 -2 -3
有没有更好的方法对此进行编码?
注意:我在论文中使用下标j代替i
您只需避免创建即可获得不错的加速效果rep(0, (kk-1)*(kk-2)/2)
。删除此步骤可以大大加快操作速度,这一事实使我认为如果不使用,则不可能更快Rcpp
create_Uj2 <- function(uj) {
q <- length(uj)
if (q == 1) return(0)
nr <- q
nc <- q*(q-1)/2
Uj <- matrix(0, nrow = nr, ncol = nc)
for (kk in 2:nr) {
Uj[kk, ((kk-1)*(kk-2)/2 + 1):(kk*(kk-1)/2)] <- uj[1:(kk-1)]
}
-Uj
}
all.equal(create_Uj2(1:400), create_Uj(1:400))
# [1] TRUE
microbenchmark(
create_Uj2(1:400),
create_Uj(1:400),
times = 10,
unit = 'relative'
)
# Unit: relative
# expr min lq mean median uq max neval
# create_Uj2(1:400) 1.000000 1.000000 1.000000 1.000000 1.000000 1.0000000 10
# create_Uj(1:400) 2.070708 1.847242 1.529658 2.028489 1.496275 0.9441935 10
本文收集自互联网,转载请注明来源。
如有侵权,请联系 [email protected] 删除。
我来说两句