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

[Deepayan Sarkar]

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.

> 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


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