Martin Maechler
2021-Apr-17 15:38 UTC
[R] Matrix::solve() with 1-d arrays -- treating "array" as "numeric"?
>>>>> 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 at 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.
Deepayan Sarkar
2021-Apr-19 04:26 UTC
[R] Matrix::solve() with 1-d arrays -- treating "array" as "numeric"?
On Sat, Apr 17, 2021 at 9:08 PM Martin Maechler <maechler at 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 at 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.