Matrix::solve() with 1-d arrays

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

Matrix::solve() with 1-d arrays

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

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

?

-Deepayan

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: 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

    > ______________________________________________
    > [hidden email] 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.

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

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

Deepayan Sarkar
On Sat, Apr 17, 2021 at 9:08 PM Martin Maechler
<[hidden email]> 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
>
>     > ______________________________________________
>     > [hidden email] 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.

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: 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
    > <[hidden email]> 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

______________________________________________
[hidden email] 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.