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

[Martin Maechler]

>>>>> 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",

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



    > -Deepayan

    > ______________________________________________
    > R-help using r-project.org mailing list -- To UNSUBSCRIBE and
    > more, see https://stat.ethz.ch/mailman/listinfo/r-help
    > PLEASE do read the posting guide
    > http://www.R-project.org/posting-guide.html and provide
    > commented, minimal, self-contained, reproducible code.