Multiply matrices to yield 3-dimension array (or, dot product but `max` instead of `sum`)

I've spent most of my R time in the tidyverse, and am now trying to branch out into more matrix math – but I'm very unfamiliar.

Perhaps it its core, my question is "is apply the 'right'/'most current'/'best' function for generically working with matrices and arrays?" Is there a standard arrays package, like dplyr for data frames? Or are there matrix- and array-oriented functions I should be aware of.

But let's get into it

Right now I have two matrices. I would like to multiply the two matrices, as in a dot product, but rather than summing the results of each vector multiplication, I want to take the max. When multiplying a 5x5 matrix by a 5x2 matrix, this should result in a 5x2 matrix.

I already have a very simple apply solution for this, but I'm not sure it's taking the most advantage of R's vectorization.

B <- matrix(c(0,1,1,0,0,
              1,0,0,1,0,
              1,0,0,0,1,
              0,1,0,0,1,
              0,0,1,1,0),
            ncol = 5, dimnames = list(letters[1:5], letters[1:5]))

M <- data.frame(v1 = c(32, 1, 0, 14, 3), 
                v2 = c( 2, 6, 4,  9, 7))
rownames(M) <- letters[1:5]
M <- as.matrix(M)

My apply solution

So far, I've done using apply to break matrix B into individual rows, and multiple each by M. I then use matrixStats::colMaxs to take the maximum values of the output.

apply(B, 1, \(i) colMaxs(i * M)) %>% t()
#   v1 v2
# a  1  6
# b 32  9
# c 32  7
# d  3  7
# e 14  9

Alternately, I could have delayed taking the column maxes, first making a 3-dimensional array:

B %>% 
  apply(1, \(i) i * M) %>%
  array(dim = c(5, 2, 5),
        dimnames = list(rownames(M), colnames(M), rownames(M))) %>%
  apply(3, colMaxs)

Is there a more proper way to create this array from matrices B and M?

I imagine there's also a version where I don't make matrices but leave everything as vectors, appropriately repeat each one, and then create an array at the end, but if so, surely there's a package to help this confusing dimension-tracking:

# This yields the wrong result
rep(as.vector(B), 2) * rep(as.vector(M), 5) %>%
  array(dim = c(5, 2, 5)) %>%
  apply(3, colMaxs)

Question

Is there a non-apply way to multiply matrices, resulting in a 3-dimensional array, where each element is the product of an element-wise multiplication? Should I be looking for one, or is apply the appropriate tool for the job?

I prefer your initial apply approach, its explicit about what it does, and it doesnt introduce explicit constructions with room for errors like you noticed can be an issue

I could be wrong but I think you may be mistaken in describing what you are doing as element_wise; it appears to be non-elementwise from my understanding; you need to have processed an entire column and row before having knowledge of a resulting cell value; rather than being able to look at say 2 cells and knowing what the resulting cell value would be.

Also you originally propose a 2 dimensional result; one of your approaches to achieving that was a concept involving starting with a 3 dimensional matrix; but at the end of your post you write as though your goal is to constuct a 3 dimensional array, in contrast to your opening. this creates ambiguity.

Thank you for your comments, and yes, I suppose that I was confusing. For this particular application I wanted a 2D final result, but I’m interested in the more generic route to get there, through the 3D array.

In my actual use case, I am trying to take the max, min, mean, and sum of each resulting product vector. I could fit those into the first apply solution, or I could run them as matrixStats functions on the resulting larger array.