[Rd] preserving sparse matrices (Matrix)
Martin Maechler
maechler at stat.math.ethz.ch
Wed Dec 20 15:43:30 CET 2006
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
More information about the R-devel
mailing list