[R] %*% in Matrix objects

[Martin Maechler]

>>>>> "Jose" == Jose Quesada <quesada at gmail.com>
>>>>>     on Sat, 27 Jan 2007 23:42:34 +0100 writes:

    Jose> Hi Martin, Thanks for your detailed answer.

    Jose> x <- Matrix(1:12, 3,4, sparse = TRUE)

    >> I hope that you are aware of the fact that it's not
    >> efficient at all to store a dense matrix (it has *no* 0
    >> entry) as a sparse one..
    >> 
    >> and your posting is indeed an incentive for the Matrix
    >> developers to improve that part ... ;-)
    >> 

    Jose> Yes, the toy example is not sparse but the actual data
    Jose> is, and very large; I'm aware that coercing a dense
    Jose> matrix into the Sparse format is not leading to any
    Jose> saving (on the contrary). I'm talking about a real
    Jose> application with large sparse matrices; from now on,
    Jose> I'll post small examples using sparse matrices as well
    Jose> to avoid confusion.

ok.

      Jose> so I tried

      Jose> x = matrix(1:12,3,4)
      Jose> x = as(x, "CsparseMatrix")
      Jose> xnorms  = sqrt(colSums(x^2))
      Jose> xnorms = as(xnorms, "CsparseMatrix")
      Jose> (xnormed = t(x) * (1/xnorms))

      Jose> But now, instead of a warning I get
      Jose> "Error: Matrices must have same dimensions in t(x) * (1/xnorms)"

    >> yes.  And the same happens with traditional matrices -- and well so:
    >> For arithmetic with matrices (traditional or "Matrices"),
    >> 
    >> A o B       (o in {"+", "*", "^", ....})
    >> -----
    >> 
    >> does require that matrices A and B are ``conformable'', i.e.,
    >> have exact same dimensions.
    >> 
    >> Only when one of A or B is *not* a matrix,
    >> then the usual S-language recycling rules are applied,
    >> and that's what you were using in your first example
    >> (<Matrix> * <numeric>) above.
    >> 

    Jose> Right. So this means that the * operator is not
    Jose> overloaded in Matrix (that is, if I use it, I'll get
    Jose> my Matrix coherced to matrix. Is that correct?

no.  The "*"  is overloaded (read on)

    Jose> Does this mean that there is no easy way to do element-by-element
    Jose> multiplication without leaving the sparse Matrix format?

No. There is an easy way:
If you multiply (or add or ..) two sparse matrices of matching dim(), the
result will be sparse. Also if use a "scalar" (length-1 vector)
with a Matrix, the result remains sparse (where appropriate) :

  > (x <- Matrix(c(0,1,0,0), 3,3))
  3 x 3 sparse Matrix of class "dtCMatrix"

  [1,] . . .
  [2,] 1 . .
  [3,] . 1 .
  Warning message:
  data length [4] is not a sub-multiple or multiple of the number of rows [3] in matrix 
  > (2 * x) + t(x)
  3 x 3 sparse Matrix of class "dgCMatrix"

  [1,] . 1 .
  [2,] 2 . 1
  [3,] . 2 .
  > ((2 * x) + t(x)) * t(x)
  3 x 3 sparse Matrix of class "dgCMatrix"

  [1,] . 1 .
  [2,] . . 1
  [3,] . . .


What you tried to do,  <sparse> * <vector-of-length-gt-1>, will
only result in a sparse matrix in the next version of the Matrix
package. 

    Jose> I suspect I'm facing the drop=T as before...
    >> why??

    Jose> Because when I got a row out of a Matrix object, the
    Jose> resulting vector is not of class Matrix but numeric,
    Jose> and then (<Matrix> * <numeric>) is applied.


    Jose> Last, I shouldn't consider myself the most standard
    Jose> user of the matrix package, since my lineal algebra is
    Jose> really basic. But in any case, you should know that
    Jose> your package is being enormously useful for me. Keep
    Jose> up the good work. And if I can help by posting my very
    Jose> basic questions, I'm glad to help.

Ok, thanks for the flowers :-)
Martin