在 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 和一些小例子。

这里还有一些选项。 maxmin1maxmin2 实现 max-min 组合。 maxprod1maxprod2 实现 max-product 组合。

maxmin1maxprod1 将与嵌套的 for 循环执行类似(如果不差的话),因为 apply 在其主体中有一个循环。

maxmin2maxprod2 是仅依赖向量化函数的优化版本。当 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