[R] Matrix::solve() with 1-d arrays -- treating "array" as "numeric"?
Deepayan Sarkar
deep@y@n@@@rk@r @end|ng |rom gm@||@com
Mon Apr 19 06:26:58 CEST 2021
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.
More information about the R-help
mailing list