What is the simplest way to implement fuzzy relational composition of two matrices in R?

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.

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.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