在 R 中实现两个矩阵的模糊关系组合的最简单方法是什么?
What is the simplest way to implement fuzzy relational composition of two matrices in R?
在R中实现两个矩阵的模糊关系组合的最简单方法是什么?我编写了它的一个版本,但据说速度很慢,所以我想知道是否有矢量化操作可以使其更快
circ_prod <- function(R,S) {
if(ncol(R) != nrow(S)) errorCondition("dimensions don't match")
R_circ_S <- matrix(0, nrow(R), ncol(S))
for (i in 1:nrow(R)) {
for (k in 1:ncol(S))
R_circ_S[i,k] <- max(pmin(R[i,],S[,k]))
}
R_circ_S
}
有一个 link 解释了为什么这样做很重要 what is fuzzy relational composition 和一些小例子。
这里还有一些选项。 maxmin1
和 maxmin2
实现 max-min 组合。 maxprod1
和 maxprod2
实现 max-product 组合。
maxmin1
和 maxprod1
将与嵌套的 for
循环执行类似(如果不差的话),因为 apply
在其主体中有一个循环。
maxmin2
和 maxprod2
是仅依赖向量化函数的优化版本。当 ncol(R)
超过 2 时,他们使用包 matrixStats
中的 colMaxs
来有效地找到列最大值。
在 C 中实现这些函数可以进一步优化。
maxmin1 <- function(R, S) {
t(apply(R, 1L, function(x) apply(pmin(S, x), 2L, max)))
}
maxmin2 <- function(R, S) {
m <- (d <- dim(R))[1L]
p <- d[2L]
n <- dim(S)[2L]
if (p == 1L) {
return(matrix(pmin.int(rep(R, each = n), S), m, n, byrow = TRUE))
}
r <- sequence.default(nvec = rep.int(p, m * n), from = rep(seq_len(m), each = n), by = m)
x <- pmin.int(S, R[r])
if (p == 2L) {
y <- x[seq.int(1L, length(x), 2L) + (x[c(TRUE, FALSE)] < x[c(FALSE, TRUE)])]
} else {
y <- matrixStats::colMaxs(matrix(x, p))
}
matrix(y, m, byrow = TRUE)
}
maxprod1 <- function(R, S) {
t(apply(R, 1L, function(x) apply(S * x, 2L, max)))
}
maxprod2 <- function(R, S) {
m <- (d <- dim(R))[1L]
p <- d[2L]
n <- dim(S)[2L]
if (p == 1L) {
return(matrix(rep(R, each = n) * S, m, n, byrow = TRUE))
}
r <- sequence.default(nvec = rep.int(p, m * n), from = rep(seq_len(m), each = n), by = m)
x <- as.double(S) * R[r]
if (p == 2L) {
y <- x[seq.int(1L, length(x), 2L) + (x[c(TRUE, FALSE)] < x[c(FALSE, TRUE)])]
} else {
y <- matrixStats::colMaxs(matrix(x, p))
}
matrix(y, m, byrow = TRUE)
}
R <- matrix(c(0.7, 0.8, 0.6, 0.3), 2L, 2L)
R
## [,1] [,2]
## [1,] 0.7 0.6
## [2,] 0.8 0.3
S <- matrix(c(0.8, 0.1, 0.5, 0.6, 0.4, 0.7), 2L, 3L)
S
## [,1] [,2] [,3]
## [1,] 0.8 0.5 0.4
## [2,] 0.1 0.6 0.7
maxmin1(R, S)
## [,1] [,2] [,3]
## [1,] 0.7 0.6 0.6
## [2,] 0.8 0.5 0.4
maxmin2(R, S)
## [,1] [,2] [,3]
## [1,] 0.7 0.6 0.6
## [2,] 0.8 0.5 0.4
microbenchmark::microbenchmark(maxmin1(R, S), maxmin2(R, S), times = 10000L)
## Unit: microseconds
## expr min lq mean median uq max neval
## maxmin1(R, S) 34.645 36.244 40.519689 36.695 37.884 3165.159 10000
## maxmin2(R, S) 3.198 3.567 4.236501 3.813 4.018 1030.412 10000
maxprod1(R, S)
## [,1] [,2] [,3]
## [1,] 0.56 0.36 0.42
## [2,] 0.64 0.40 0.32
maxprod2(R, S)
## [,1] [,2] [,3]
## [1,] 0.56 0.36 0.42
## [2,] 0.64 0.40 0.32
microbenchmark::microbenchmark(maxprod1(R, S), maxprod2(R, S), times = 10000L)
## Unit: microseconds
## expr min lq mean median uq max neval
## maxprod1(R, S) 25.543 26.814 29.412076 27.265 28.003 1174.527 10000
## maxprod2(R, S) 2.911 3.239 3.723411 3.403 3.608 932.586 10000
在R中实现两个矩阵的模糊关系组合的最简单方法是什么?我编写了它的一个版本,但据说速度很慢,所以我想知道是否有矢量化操作可以使其更快
circ_prod <- function(R,S) {
if(ncol(R) != nrow(S)) errorCondition("dimensions don't match")
R_circ_S <- matrix(0, nrow(R), ncol(S))
for (i in 1:nrow(R)) {
for (k in 1:ncol(S))
R_circ_S[i,k] <- max(pmin(R[i,],S[,k]))
}
R_circ_S
}
有一个 link 解释了为什么这样做很重要 what is fuzzy relational composition 和一些小例子。
这里还有一些选项。 maxmin1
和 maxmin2
实现 max-min 组合。 maxprod1
和 maxprod2
实现 max-product 组合。
maxmin1
和 maxprod1
将与嵌套的 for
循环执行类似(如果不差的话),因为 apply
在其主体中有一个循环。
maxmin2
和 maxprod2
是仅依赖向量化函数的优化版本。当 ncol(R)
超过 2 时,他们使用包 matrixStats
中的 colMaxs
来有效地找到列最大值。
在 C 中实现这些函数可以进一步优化。
maxmin1 <- function(R, S) {
t(apply(R, 1L, function(x) apply(pmin(S, x), 2L, max)))
}
maxmin2 <- function(R, S) {
m <- (d <- dim(R))[1L]
p <- d[2L]
n <- dim(S)[2L]
if (p == 1L) {
return(matrix(pmin.int(rep(R, each = n), S), m, n, byrow = TRUE))
}
r <- sequence.default(nvec = rep.int(p, m * n), from = rep(seq_len(m), each = n), by = m)
x <- pmin.int(S, R[r])
if (p == 2L) {
y <- x[seq.int(1L, length(x), 2L) + (x[c(TRUE, FALSE)] < x[c(FALSE, TRUE)])]
} else {
y <- matrixStats::colMaxs(matrix(x, p))
}
matrix(y, m, byrow = TRUE)
}
maxprod1 <- function(R, S) {
t(apply(R, 1L, function(x) apply(S * x, 2L, max)))
}
maxprod2 <- function(R, S) {
m <- (d <- dim(R))[1L]
p <- d[2L]
n <- dim(S)[2L]
if (p == 1L) {
return(matrix(rep(R, each = n) * S, m, n, byrow = TRUE))
}
r <- sequence.default(nvec = rep.int(p, m * n), from = rep(seq_len(m), each = n), by = m)
x <- as.double(S) * R[r]
if (p == 2L) {
y <- x[seq.int(1L, length(x), 2L) + (x[c(TRUE, FALSE)] < x[c(FALSE, TRUE)])]
} else {
y <- matrixStats::colMaxs(matrix(x, p))
}
matrix(y, m, byrow = TRUE)
}
R <- matrix(c(0.7, 0.8, 0.6, 0.3), 2L, 2L)
R
## [,1] [,2]
## [1,] 0.7 0.6
## [2,] 0.8 0.3
S <- matrix(c(0.8, 0.1, 0.5, 0.6, 0.4, 0.7), 2L, 3L)
S
## [,1] [,2] [,3]
## [1,] 0.8 0.5 0.4
## [2,] 0.1 0.6 0.7
maxmin1(R, S)
## [,1] [,2] [,3]
## [1,] 0.7 0.6 0.6
## [2,] 0.8 0.5 0.4
maxmin2(R, S)
## [,1] [,2] [,3]
## [1,] 0.7 0.6 0.6
## [2,] 0.8 0.5 0.4
microbenchmark::microbenchmark(maxmin1(R, S), maxmin2(R, S), times = 10000L)
## Unit: microseconds
## expr min lq mean median uq max neval
## maxmin1(R, S) 34.645 36.244 40.519689 36.695 37.884 3165.159 10000
## maxmin2(R, S) 3.198 3.567 4.236501 3.813 4.018 1030.412 10000
maxprod1(R, S)
## [,1] [,2] [,3]
## [1,] 0.56 0.36 0.42
## [2,] 0.64 0.40 0.32
maxprod2(R, S)
## [,1] [,2] [,3]
## [1,] 0.56 0.36 0.42
## [2,] 0.64 0.40 0.32
microbenchmark::microbenchmark(maxprod1(R, S), maxprod2(R, S), times = 10000L)
## Unit: microseconds
## expr min lq mean median uq max neval
## maxprod1(R, S) 25.543 26.814 29.412076 27.265 28.003 1174.527 10000
## maxprod2(R, S) 2.911 3.239 3.723411 3.403 3.608 932.586 10000