Re: [Rd] preserving sparse matrices (Matrix)

From: Martin Maechler <maechler_at_stat.math.ethz.ch>
Date: Wed 20 Dec 2006 - 14:43:30 GMT

Thank you, Tamas,

>>>>> "Tamas" == Tamas K Papp <tpapp@princeton.edu> >>>>> on Tue, 19 Dec 2006 18:25:16 +0100 writes:

    Tamas> Hi,
    Tamas> I have sparse (tridiagonal) matrices, and I use the Matrix package for
    Tamas> handling them.  For certain computations, I need to either set the
    Tamas> first row to zero, or double the last row.  I find that neither
    Tamas> operation preserves sparsity, eg

>> m <- Diagonal(5)
>> m

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

showClass("ddiMatrix")

will show you that 'm' is not a 'sparseMatrix' but a 'diagonalMatrix' which inherits from 'denseMatrix'; since I had decided that 'sparseMatrix' should be confined to those matrices that store "location + content"; This is in line with the fact that a packed triangular matrix also inherits from 'denseMatrix'.
(Also, to store a diagonalMatrix as sparseMatrix leads to larger  object size).

I agree however that this is not obvious to the casual user and also that you should not *have* to know this.

>> m[1,] <- 0
>> m

    Tamas> 5 x 5 Matrix of class "dgeMatrix"
    Tamas> [,1] [,2] [,3] [,4] [,5]
    Tamas> [1,]    0    0    0    0    0
    Tamas> [2,]    0    1    0    0    0
    Tamas> [3,]    0    0    1    0    0
    Tamas> [4,]    0    0    0    1    0
    Tamas> [5,]    0    0    0    0    1

>> m <- Diagonal(5)
>> rbind(m,1:5)
    Tamas> 6 x 5 Matrix of class "dgeMatrix"
    Tamas> [,1] [,2] [,3] [,4] [,5]
    Tamas> [1,]    1    0    0    0    0
    Tamas> [2,]    0    1    0    0    0
    Tamas> [3,]    0    0    1    0    0
    Tamas> [4,]    0    0    0    1    0
    Tamas> [5,]    0    0    0    0    1
    Tamas> [6,]    1    2    3    4    5

Yes, this is an infelicity (and I'd say for the first example even a bug) in the 'Matrix' package which I will fix for the next version.

    Tamas> Currently, I am calling Matrix to make these sparse again, but my
    Tamas> matrices are large (100x100).  Is there a way to perform the
    Tamas> operations I want without losing sparsity?

Yes --- however, not as easily as I thought for the current
	version of 'Matrix' -- due to more infelicities... :

Here's a workable script {which has output as comments}:

m0 <- Diagonal(5)

### 1st problem : ----------------------------

rbind(m0, 1:5)# -> dense

## Need this for now {not anymore in Matrix 0.9975-8 } : .rb2m <- function(x, y) callGeneric(x, matrix(y, nrow = 1))

setMethod("rbind2", signature(x = "Matrix",    y = "numeric"), .rb2m)
setMethod("rbind2", signature(x = "ddiMatrix", y = "numeric"), .rb2m)
setMethod("rbind2", signature(x = "ddiMatrix", y = "matrix"),
          function(x,y) rbind2(suppressWarnings(as(x,"CsparseMatrix")),
                               Matrix:::as_Csparse(y)))

rbind(m0, 1:5)
## 6 x 5 sparse Matrix of class "dgCMatrix"

## [1,] 1 . . . .
## [2,] . 1 . . .
## [3,] . . 1 . .
## [4,] . . . 1 .
## [5,] . . . . 1
## [6,] 1 2 3 4 5

### 2nd problem: ----------------------------

diag2Sparse <- function(x) ## this is a workaround for a current "Matrix" bug

    Matrix:::diagU2N(suppressWarnings(as(x, "CsparseMatrix")))

m1 <- diag2Sparse(m0)

m1[1,] <- 0 ; m1
m1[,3] <- 3 ; m1
m1[1:3, 3] <- 0

m1
Matrix:::drop0(m1)# remove the "0"
## Should drop0() be exported and documented? -- probably yes!

I hope this helps.

Next time, we (Prof Douglas Bates and me) will appreciate if you approach us directly.

Best regards,
Martin



R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel Received on Thu Dec 21 12:32:02 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Thu 21 Dec 2006 - 02:31:04 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel. Please read the posting guide before posting to the list.