1

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.

monotonic
  • 394
  • 4
  • 20
  • Sounds interesting, but we need to see sample input, expected output, and some kind of explanation of how to get from one to the other. – Mikael Jagan Feb 17 '22 at 03:57
  • Should be in the link @MikaelJagan – monotonic Feb 17 '22 at 03:58
  • 3
    Hi @monotonic. This is not a good question for SO in its current form. SO is about *specific* code issues, that can be demonstrated & verified with sample data and some attempt at translating the problem into code. At the moment your question either reads like you're asking somebody else to do a project for you, or you're asking for a package /tool which is OT here on SO. I suggest editing your post addressing Mikael Jagan's suggestions and then flag the post for "Reopen". – Maurits Evers Feb 17 '22 at 04:08
  • 1
    Consider that, if your link breaks, then any answers will be entirely without context for someone who comes across this question in the future. – Mikael Jagan Feb 17 '22 at 04:27
  • @MikaelJagan that's a good point and thanks for mentioning that. I just coded a for loop version, and thanks for providing the answer. Can we reopen this questions so others may see it too? – monotonic Feb 17 '22 at 05:02
  • @MauritsEvers good to reopen? – monotonic Feb 17 '22 at 05:06
  • @monotonic Flagged for reopen. Needs one more vote. – Maurits Evers Feb 17 '22 at 05:38

1 Answers1

2

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
Mikael Jagan
  • 9,012
  • 2
  • 17
  • 48