R help - Jan 2007 - %*% in Matrix objects

Dear R users,

I need to normalize a bunch of row vectors. At a certain point I need to divide
a matrix by a vector of norms. I find that the behavior of Matrix objects
differs from normal matrix objects. Example the following code examples differ
only in xnormed changing from normal to Matrix object:

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

This produces a "warning: coercing sparse to dense matrix for arithmetic
in: t(x) * 1/xnorms." but gets the result (a 4 x 3 matrix)

I want to stay in sparse format anyway (if it helps!) so I tried

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

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

If I transpose the norms, the error dissapears, but the result is 1 x 4 (not 3 x
4 as before).

I suspect I'm facing the drop=T as before...
Also, it seems that in normal matrix objects %*% behaves the same as *, but in
Matrix objects that is not the case.

What am I missing?

-- 
Thanks,
-Jose

--
Jose Quesada, PhD
Research fellow, Psychology Dept.
Sussex University, Brighton, UK
http://www.andrew.cmu.edu/~jquesada

>>>>> "Jose" == Jose Quesada <quesada at
gmail.com>
>>>>>     on Fri, 26 Jan 2007 05:24:12 +0100 writes:
    Jose> Dear R users,
    Jose> I need to normalize a bunch of row vectors. At a certain point I
need to divide a matrix by a vector of norms. I find that the behavior of Matrix
objects differs from normal matrix objects.

I believe you are showing evidence for that; though I know it's
still true, and will be less true for the next release of the
Matrix package.

    Jose> Example the following code examples differ only in
    Jose> xnormed changing from normal to Matrix object:

    Jose> x = matrix(1:12,3,4)
    Jose> x = as(x, "CsparseMatrix")

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

    Jose> xnorms  = sqrt(colSums(x^2))
    Jose> (xnormed = t(x) * (1/xnorms))

    Jose> This produces a "warning: coercing sparse to dense matrix for
arithmetic
    Jose> in: t(x) * 1/xnorms." but gets the result (a 4 x 3 matrix)

    Jose> I want to stay in sparse format anyway
    Jose>  (if it helps!) 

what should it help for? Are you talking about a real
application with a proper sparse matrix as opposed to the toy
example here? In that case I agree, and as a matter of fact, the
source code of the Matrix package leading to the above warning
is the following

	      } else {
		  ## FIXME: maybe far from optimal:
		  warning("coercing sparse to dense matrix for arithmetic")
		  callGeneric(as(e1, "dgeMatrix"), e2)
	      }

and your posting is indeed an incentive for the Matrix developers
to improve that part ... ;-)

    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> If I transpose the norms, the error dissapears, but the result is 1
x 4 (not 3 x 4 as before).

That would be a bug of *not* giving an error... but I can't see
it. Can you please give an exact example -- as you well gave otherwise; and 
BTW, do you use a sensible sparse matrix such as
 (x <- Matrix(c(0,0,1,0), 3,4))


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

    Jose> Also, it seems that in normal matrix objects %*%
    Jose> behaves the same as *,

not at all!!  If that was the case, the inventors of the S
language would never have introduced '%*%' !

    Jose> but in Matrix objects that is not the case.

    Jose> What am I missing?

I guess several things, see above, notably how basic
matrix+vector - arithmetic is defined to in the S language.

Regards,
Martin Maechler, ETH Zurich

Hi Martin,
Thanks for your detailed answer.

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

	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.
>
Right. So this means that the * operator is not overloaded in Matrix (that is,
if I use it, I'll get my Matrix coherced to matrix. Is that correct?
Does this mean that there is no easy way to do element-by-element
multiplication
without leaving the sparse Matrix format?

    Jose> I suspect I'm facing the drop=T as
before...
> why??
Because when I got a row out of a Matrix object, the resulting vector is not of
class Matrix but numeric, and then (<Matrix> * <numeric>) is
applied.

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



-- 
Cheers,
-Jose

--
Jose Quesada, PhD
Research fellow, Psychology Dept.
Sussex University, Brighton, UK
http://www.andrew.cmu.edu/~jquesada