[Rd] preserving sparse matrices (Matrix)

[Martin Maechler]

Thank you, Tamas,

>>>>> "Tamas" == Tamas K Papp <tpapp at 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