[R] Matrix::solve() with 1-d arrays -- treating "array" as "numeric"?

[Martin Maechler]

>>>>> Deepayan Sarkar 
>>>>>     on Mon, 19 Apr 2021 09:56:58 +0530 writes:

    > On Sat, Apr 17, 2021 at 9:08 PM Martin Maechler
    > <maechler using stat.math.ethz.ch> wrote:
    >> 
    >> >>>>> Deepayan Sarkar >>>>> on Fri, 16 Apr 2021 11:34:20
    >> +0530 writes:
    >> 
    >> > I get what I initially thought was unexpected behaviour
    >> from:
    >> 
    >> > x <- tapply(runif(100), sample(5, 100, TRUE), mean) >
    >> solve(Diagonal(5), x) > # Error: not-yet-implemented
    >> method for solve(<ddiMatrix>, <array>).  > # ->> Ask the
    >> package authors to implement the missing feature.
    >> 
    >> ((why did you not ask the package authors ?? --- never
    >> mind))
    >> 
    >> 
    >> > This is because x is a 1-D array, so the operation is
    >> not > well-defined. Would it make sense for Matrix to
    >> support this (treat > 1-D arrays as column vectors, as it
    >> does for plain vectors)? Or should > I make my intent
    >> clear with
    >> 
    >> > solve(Diagonal(5), as.vector(x))
    >> 
    >> well ...
    >> 
    >> The "fun" thing is that it actually works when Matrix
    >> methods are not fully available, i.e., if you do *not* do
    >> require(Matrix) or equivalent, but rather only load the
    >> Matrix namespace via
    >> 
    >> solve(Matrix::Diagonal(5), x)
    >> 
    >> actually currently works correctly by some "good
    >> coincidence" (I have not yet tried to understand, as
    >> that's a bit painful: selectMethod("solve",
    >> c("ddiMatrix", "array")) is "lying" here ! )
    >> 
    >> However this looks like a more general problem with S4
    >> methods -- and probably a good reason for asking on
    >> R-help -- namely, the fact that d1-dimensional (numeric)
    >> arrays are not automatically treated as (numeric) vectors
    >> i.e. class "numeric" wrt S4 methods.
    >> 
    >> In the following case the solve() - coincidence does not
    >> help, BTW.
    >> 
    >> Diagonal(3) %*% array(1:3)
    >> 
    >> ## Error in Diagonal(3) %*% array(1:3) : ##
    >> not-yet-implemented method for <ddiMatrix> %*% <array>
    >> 
    >> 
    >> In principle, we should consider a way to tell that
    >> "array" should be tried as "vector",

    > Actually, if you think compatible 1-d numeric arrays
    > should 'work', then all I was thinking of was something
    > along the following lines: Add an additional

    > setMethod("solve", c("ANY", "array"), function(a, b, ...)
    > ...

    > which would basically do a dimension check for b, for 1-d
    > numeric arrays call solve(a, as.vector(b), ...), and error
    > out for dim(b) > 2.  The actual details may be more
    > involved, but that's the basic idea.

Well, of course, it's just that a method like the one above
would have to be provided for all signatures where something
corresponding for "numeric" now exists ..
and these may easily get to many dozens at least.

In other words  solve() is just one special case of very many
and I was rather interested in finding a method to solve "all
cases" at once - "automagically".

Otherwise, supporting  1d "array" to work Matrix-pkg matrices
would not only mean providing these many dozen of methods now,
but
also for every new functionality that uses S4 methods with
numeric vectors to have to always provide an additional "array"
method, i.e.,  increasing the Matrix-maintenance burden yet more.




    >> possibly via something like setIs("array", "vector") or
    >> rather setIs("array", "numeric") because in the Matrix
    >> package the vectors encountered are really numeric
    >> vectors.
    >> 
    >> .. OTOH, in all of the 3 packages I co-author and which
    >> use S4 heavily, Matrix, Rmpfr, lme4, I had till now
    >> decided *not* to setIs() because it never worked as I
    >> intended, or rather had unpleasant side effects.
    >> 
    >> Here, setIs("array", "numeric", test=is.numeric)
    >> 
    >> gives
    >> 
    >> Error in setIs("array", "numeric", test = is.numeric) :
    >> cannot create a 'setIs' relation when neither of the
    >> classes (“array” and “numeric”) is local and modifiable
    >> in this package
    >> 
    >> A more successful alternative had been to use
    >> setClassUnion(), so I could consider defining
    >> 
    >> setClassUnion("mVector", c("numeric", "array"))
    >> 
    >> and replace "numeric" in many of the method signatures by
    >> "mVector" (but that would then also dispatch for 1d
    >> character arrays ... not so nicely).

    > But you already have that problem, I think:

    >> s = matrix(letters[1:10], 5, 2) solve(Diagonal(5), s)
    > Error in .M.kind(data) : not yet implemented for matrix
    > with typeof character

    > whereas

    > m = matrix(1:10, 5, 2)

    > works nicely. Unfortunately, both m and s have the same
    > class (c("matrix", "array")), so I don't think method
    > dispatch would be able to distinguish between them with
    > the current design, and you anyway need to check in the
    > solve method for c("diagonalMatrix", "matrix").

    > Best, -Deepayan

That's a convincing argument, I concur,
thank you.

So maybe such a "global" class union containing "array", and
replacing "numeric" in Matrix' current setMethod()s' signatures
would solve the problem  quite nicely (well, as "nicely" as it
can easily get).

Something we should really try (as soon as current R-forge
version of Matrix, 1.3-3, will have been on CRAN).

Note BTW, that since R 4.0.x, we have  the function

 .class2()  which provides the full
 	    "method-dispatch class vector", 
e.g.,

> .class2(array(1:3))
[1] "array"   "integer" "numeric"
> .class2(matrix(letters, 2))
[1] "matrix"    "array"     "character"
>

Martin