How can I write a `%*%` method for a base matrix S3 subclass?

Question

I'd like to write a %*% method for a base matrix subclass. My subclass is an S3 class and the documentation of help("%*%") says that %*% is a S4 generic and that S4 methods need to be written for a function of two arguments named x and y. I've written methods for S3 classes for S4 generic methods before using methods::setOldClass() and methods::setMethod() and I've looked at the source code for {Matrix} package for inspiration but for some reason I can't quite get this to work in this case. methods::showMethods() shows my method exists for my target signature but my method never seems to actually get called by R when I try to use it.

x <- diag(3)
class(x) <- c("atm2d", class(matrix()))
print(x)

     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1
attr(,"class")
[1] "atm2d"  "matrix" "array" 

The default %*% gets rid of my class attribute and I would like to preserve it.

print(x %*% x)

     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1

I've tried to create a %*% method for my class which doesn't get rid of my class attribute:

as.matrix.atm2d <- function(x, ...) {
    class(x) <- NULL
    x
}
matmult <- function(x, y) {
    v <- as.matrix(x) %*% as.matrix(y)
    class(v) <- c("atm2d", class(matrix()))
    v
}
methods::setOldClass("atm2d")
# this alternate `setOldClass()` also doesn't work
# methods::setOldClass(c("atm2d", class(matrix()))) 
methods::setMethod("%*%", 
                   c(x = "atm2d", y = "atm2d"), 
                   function(x, y) matmult(x, y))

showMethods() seems to show an S4 method with the expected signature was created:

showMethods("%*%", class = "atm2d")

Function: %*% (package base)
x="atm2d", y="atm2d"

However this method doesn't actually seem to be called by %*%:

print(x %*% x)

     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1

If my method had been called I would have instead expected it to also print out its class:

print(matmult(x, x))

     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1
attr(,"class")
[1] "atm2d"  "matrix" "array"

Answer 1

The %*% operator is internally generic, which means that dispatch happens in C code. Currently (i.e., in R 4.2.3), the corresponding C function do_matprod (defined here) contains this check:

if (PRIMVAL(op) == 0 && /* %*% is primitive, the others are .Internal() */
   (IS_S4_OBJECT(x) || IS_S4_OBJECT(y))
   && R_has_methods(op)) {
    SEXP s, value;
    /* Remove argument names to ensure positional matching */
    for(s = args; s != R_NilValue; s = CDR(s)) SET_TAG(s, R_NilValue);
    value = R_possible_dispatch(call, op, args, rho, FALSE);
    if (value) return value;
}

If neither x nor y is an S4 object, as in your example, then do_matprod proceeds to handle them as traditional matrices, not looking at the class attribute of either argument. The section of help("%*%") to which you refer:

This operator is S4 generic but not S3 generic. S4 methods need to be written for a function of two arguments named x and y.

tries to express that, but it is not particularly clear. (You did define an S4 method, after all.)

There are two main issues here:

• setOldClass allows you to define S4 methods with S3 classes in the signature, but internally generic functions only look for S4 methods when one of the arguments is an S4 object (for speed).

• %*% is not S3 generic, so even if you registered an S3 method like %*%.zzz, it would not be dispatched.

Having said that, the %*% operator will become S3 generic starting in R 4.3.0, due to be released on April 21. You'll find this entry in the NEWS:

The matrix multiply operator %*% is now an S3 generic, belonging to new group generic matrixOps. From Tomasz Kalinowski's contribution in PR#18483.

When that happens, %*% will behave like + and other internally generic members of the Ops group, in that S3 methods %*%.zzz will be dispatched where appropriate. However, S4 methods will still not be dispatched when neither argument is an S4 object.

Arguably, setMethod should be changed to signal a warning or error when you try to define S4 methods that will never be dispatched, as you have in your example.

---

Appendix

It might help to enumerate the types of generic functions and the types of methods that they dispatch, restricting attention to S3 and S4. We'll use this script to define objects for our tests, which should each be run in a new R process:

## objects.R
w <- structure(0, class = "a")
x <- structure(0, class = "b")

setOldClass("c")
y <- structure(0, class = "c")

setClass("d", contains = "numeric")
z <- new("d", 0)

Non-internally S3 generic functions

These dispatch S3 methods for S3 classes and S4 classes via UseMethod. When no method is found, they dispatch the default method *.default or (if that is not found) throw an error. They never dispatch S4 methods.

source("objects.R")

h <- .__h__. <- function(x) UseMethod("h")
.S3method("h", "default", function(x) "default")

.S3method("h", "a", function(x) "a3")
setMethod("h", "c", function(x) "c4")
setMethod("h", "d", function(x) "d4")

h <- .__h__. # needed to undo side effect of 'setMethod'
h
## function(x) UseMethod("h")

h(w)
## [1] "a3"
h(x)
## [1] "default"
h(y)
## [1] "default"
h(z)
## [1] "default"

Non-internally S4 generic functions

These dispatch S4 methods for S3 classes (formally defined with setOldClass) and S4 classes via standardGeneric. When no method is found, they dispatch the default method *@default. If the default method is S3 generic, then dispatch happens again, this time to any available S3 methods. However, often the default method just calls stop to signal an error.

source("objects.R")

h <- function(x) UseMethod("h")
.S3method("h", "default", function(x) "default")

.S3method("h", "c", function(x) "c3")
setMethod("h", "d", function(x) "d4")

h
## standardGeneric for "h" defined from package ".GlobalEnv"
## 
## function (x) 
## standardGeneric("h")
## <environment: 0x1044650b0>
## Methods may be defined for arguments: x
## Use  showMethods(h)  for currently available ones.

h@default
## Method Definition (Class "derivedDefaultMethod"):
## 
## function (x) 
## UseMethod("h")
## 
## Signatures:
##         x    
## target  "ANY"
## defined "ANY"

h(w)
## [1] "default"
h(x)
## [1] "default"
h(y)
## [1] "c3"
h(z)
## [1] "d4"

Internally generic functions

These are all defined in base as primitive functions. You would refer to the help pages or the source code to determine whether they are only S3 generic, only S4 generic, or both S3 and S4 generic. In the third case, S3 dispatch happens only if no suitable S4 methods are found. And, as I've explained already, S4 dispatch happens only if one of the arguments in the signature is an S4 object.

Let's take + and %*% as examples. Both are S4 generic but only + is S3 generic.

source("objects.R")

.S3method("+", "default", function(e1, e2) "default")

.S3method("+", "a", function(e1, e2) "a3")
.S3method("+", "b", function(e1, e2) "b3")
.S3method("+", "c", function(e1, e2) "c3")
.S3method("+", "d", function(e1, e2) "d3")

setMethod("+", c("c", "c"), function(e1, e2) "c4")
setMethod("+", c("d", "d"), function(e1, e2) "d4")

w + w
## [1] "a3"
x + x
## [1] "b3"
y + y
## [1] "c3"
z + z
## [1] "d4"

Here, S3 methods are dispatched. The first three results are obtained via S3 dispatch. The last result is obtained via S4 dispatch.

source("objects.R")

.S3method("%*%", "default", function(x, y) "default")

.S3method("%*%", "a", function(x, y) "a3")
.S3method("%*%", "b", function(x, y) "b3")
.S3method("%*%", "c", function(x, y) "c3")
.S3method("%*%", "d", function(x, y) "d3")

setMethod("%*%", c("c", "c"), function(x, y) "c4")
setMethod("%*%", c("d", "d"), function(x, y) "d4")

w %*% w
##      [,1]
## [1,]    0
x %*% x
##      [,1]
## [1,]    0
y %*% y
##      [,1]
## [1,]    0
z %*% z
## [1] "d4"

Here, S3 methods are not dispatched. The first three results are obtained via the internally defined default method. The last result is obtained via S4 dispatch.