What is the simplest way to implement fuzzy relational composition of two matrices in R? I coded a version of it but it's supposedly very slow, so I wonder if there's vectorized operations that can make it faster

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
}

There is a link that explains why doing so is important what is fuzzy relational composition and some small examples.

CodePudding user response：

Here are some more options. maxmin1 and maxmin2 implement the max-min composition. maxprod1 and maxprod2 implement the max-product composition.

maxmin1 and maxprod1 are going to perform similarly to (if not worse than) a nested for loop, since apply has a loop in its body.

maxmin2 and maxprod2 are optimized versions that rely solely on vectorized functions. When ncol(R) exceeds 2, they use colMaxs from package matrixStats to find columnwise maxima efficiently.

Implementing these functions in C would allow you optimize further.

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(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(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, 2L))
    }
    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, 2L))
    }
    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)
##      [,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 = 1000L)
## Unit: microseconds
##           expr    min      lq     mean  median     uq      max neval
##  maxmin1(R, S) 35.342 38.9090 50.00364 43.9725 54.694 1861.769  1000
##  maxmin2(R, S)  7.052  8.1385 10.22175  9.7990 11.521   78.064  1000

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 = 1000L)
## Unit: microseconds
##            expr    min     lq      mean median     uq      max neval
##  maxprod1(R, S) 26.937 28.454 33.684411 30.832 32.677 1664.805  1000
##  maxprod2(R, S)  3.075  3.526  3.935959  3.772  4.100   54.079  1000