[Rd] RFC: Matrix package: Matrix products (%*%, crossprod, tcrossprod) involving "nsparseMatrix" aka sparse pattern matrices

[Martin Maechler]

>>>>> Trevor Hastie <hastie at stanford.edu>
>>>>>     on Thu, 19 Mar 2015 16:03:38 -0700 writes:

    > Hi Martin
    > I got stung by this last week.
    > glmnet produces a coefficient matrix of class “dgCMatrix”
    > If a predictor matrix was created using sparseMatrix as follows,
    > one gets unexpected results, as this simple example shows.
    > My fix was easy (I always convert the predictor matrix to class “dgCMatrix” now)

    > Trevor

    >> y=Matrix(diag(4))


Considerably faster  (for larger n):   

	  Diagonal(4)

if you want a sparse matrix directly, there are

    .sparseDiagonal() 
and .symDiagonal()  
function


    >> y
    > 4 x 4 diagonal matrix of class "ddiMatrix"
    > [,1] [,2] [,3] [,4]
    > [1,]    1    .    .    .
    > [2,]    .    1    .    .
    > [3,]    .    .    1    .
    > [4,]    .    .    .    1

there's no problem with 'y' which is a "diagonalMatrix" and only
needs  O(n) storage  rather than  diag(n),
right ?

    >> z=sparseMatrix(1:4,1:4)
    >> z
    > 4 x 4 sparse Matrix of class "ngCMatrix"
            
    > [1,] | . . .
    > [2,] . | . .
    > [3,] . . | .
    > [4,] . . . |
    >> beta=as(Matrix(1:4),"dgCMatrix")
    >> y%*%beta
    > 4 x 1 sparse Matrix of class "dgCMatrix"
      
    > [1,] 1
    > [2,] 2
    > [3,] 3
    > [4,] 4
    >> z%*%beta
    > 4 x 1 sparse Matrix of class "ngCMatrix"
      
    > [1,] |
    > [2,] |
    > [3,] |
    > [4,] |
    >> 
Yes, the last one is what I consieder bogous.

Thank you, Trevor, for the feedback!
Martin


    >> On Mar 19, 2015, at 3:02 PM, Martin Maechler <maechler at stat.math.ethz.ch> wrote:
    >> 
    >> This is a Request For Comment, also BCCed to 390 package maintainers
    >> of reverse dependencies of the Matrix package.
    >> 
    >> Most users and package authors working with our 'Matrix' package will
    >> be using it for numerical computations, and so will be using
    >> "dMatrix" (d : double precision) matrix objects  M,   and indirectly, e.g., for
    >> M >= c  will also use "lMatrix" (l: logical i.e.  TRUE/FALSE/NA).
    >> All the following is  **not** affecting those numerical / logical
    >> computations.
    >> 
    >> A few others will know that we also have "pattern" matrices (purely
    >> binary: TRUE/FALSE, no NA) notably sparse ones, those "ngCMatrix" etc,
    >> all starting with "n" (from ``patter[n]``) which do play a prominent
    >> role in the internal sparse matrix algorithms, notably of the
    >> (underlying C code) CHOLMOD library in the so-called "symbolic"
    >> cholesky decomposition and other such operations. Another reason you
    >> may use them because they are equivalent to incidence matrices of
    >> unweighted (directed or undirected) graphs.
    >> 
    >> Now, as the subject says, I'm bringing up the topic of what should
    >> happen when these matrices appear in matrix multiplications.
    >> Somewhat by design, but also partly by coincidence,  the *sparse*
    >> pattern matrices multiplication in the Matrix package mostly builds on
    >> the CHOLMOD library `cholmod_ssmult()` function which implements
    >> "Boolean arithmetic" for them, instead of regular arithmetic:
    >> "+" is logical "or"
    >> "*" is  logical "and".
    >> Once we map  TRUE <-> 1  and  FALSE <-> 0, the only difference between
    >> boolean and regular arithmetic is that "1+1 = 1" in the (mapped)
    >> boolean arithmetic, because  "TRUE | TRUE" is TRUE in original logic.
    >> 
    >> The drawback of using the boolean arithmetic here is the "clash" with
    >> the usual numeric arithmetic, and arithmetic in R where logical is
    >> coerced to integer (and that to "double") when certain numerical
    >> functions/operations are used.
    >> 
    >> A more severe problem --- which I had not been aware of until
    >> relatively recently -- is the fact that  the CHOLMD function
    >> cholmod_ssdmult(A, B)
    >> treats *both* A and B as "pattern" as soon as one of them is a
    >> (sparse) pattern matrix.
    >> And this is - I say - in clear contrast to what R users would expect:
    >> If you multiply a numeric with a "kind of logical" matrix (a pattern
    >> one), you will expect that the
    >> TRUE/FALSE matrix will be treated as a 1/0 matrix because it is
    >> combined with a numeric matrix.
    >> So we could say that in this case, the Matrix package behavior is
    >> clearly bugous .... but still it has been the behavior for the last 10
    >> years or so.
    >> 
    >> RFC 1: "Change 1":
    >> I currently propose to change this behavior for the upcoming release
    >> of Matrix (version 1.2-0),  though I have no idea if dependent
    >> packages would partly fail their checks or otherwise have changed
    >> behavior subsequently.
    >> The change seems sensible, since I think if your package relied on
    >> this behavior, it was inadvertent and accidental.
    >> Still you may differ in your opinion about this change nr.1
    >> 
    >> RFC 2: "Change 2":
    >> This change would be more radical, and something I would not plan for
    >> the upcoming release of Matrix, but possibly for an update say one or
    >> two months later or so:  It concerns the matrix products when *both*
    >> matrices are pattern.  A situation where the boolean arithmetic may
    >> really make sense and where indeed packages may have depended on the
    >> current behavior  ("T + T  |--> T"). ... although that is currently
    >> only used for *sparse* pattern matrices, not for dense ones.
    >> 
    >> Further, it may still seem surprising that matrix multiplication does
    >> not behave numerically for a pair of such matrices, and by the
    >> principle of "least surprise" we should provide the boolean arithmetic
    >> matrix products in another way than  by the   standard  %*%,
    >> crossprod()  and  tcrossprod() functions.
    >> So one possibility could be to change the standard functions to behave
    >> numerically,
    >> and e.g., use   %&%  (replace the numeric "*" by a logical "&")  and
    >> crossprod(A,B, boolean=TRUE),  tcrossprod(A,B, boolean=TRUE)
    >> for the three  boolean arithmetic  version of matrix multiplications.
    >> 
    >> What do you think about this?   I'm particularly interested to hear
    >> from authors and users of  packages such as 'arules'  which IIRC
    >> explicitly work with sparse pattern matrices.
    >> 
    >> Thank you for your thoughts and creative ideas,
    >> Martin Maechler, ETH Zurich

    > ----------------------------------------------------------------------------------------
    > Trevor Hastie                                   hastie at stanford.edu <mailto:hastie at stanford.edu>  
    > Professor, Department of Statistics, Stanford University
    > Phone: (650) 725-2231                 Fax: (650) 725-8977  
    > URL: http://www.stanford.edu/~hastie <http://www-stat.stanford.edu/~hastie>  
    > address: room 104, Department of Statistics, Sequoia Hall
    > 390 Serra Mall, Stanford University, CA 94305-4065  
    > --------------------------------------------------------------------------------------