Sparse matrix no longer sparse (Matrix Package)

Thank you, John. I originally used dput() but the output is huge. However, here is a reproducible example of what I think very unexpected behavior of some matrix functions.
### Create a symmetric matrix of class dsCMatrix
A <- as(diag(5, 10), 'dsCMatrix')
A[1, 5] <- A[5,1] <- 2

### Create a diagonal matrix of class ddi
D <- Diagonal(10, 1)

### This works as it should
aa <- Cholesky(A %*% D)

### Now, let's only change one element of D to be non-integer
D[1] <- 1.5

### AD is still symmetric, but here the Cholesky function complains that it is not
aa <- Cholesky(A %*% D)
Error in Cholesky(as(A, "symmetricMatrix"), perm = perm, LDL = LDL, super = super,  : 
  error in evaluating the argument 'A' in selecting a method for function 'Cholesky': Error in asMethod(object) : 
  not a symmetric matrix; consider forceSymmetric() or symmpart()

### For fun try this
L <- update(aa, as(A %*% D, 'symmetricMatrix'))
Error in asMethod(object) : 
  not a symmetric matrix; consider forceSymmetric() or symmpart()

### This does indeed work, but should I need to implement this step?

Cholesky(forceSymmetric(A %*% D))

So, there is something about changing the elements of a ddi matrix that causes subsequent problems. Is there a good reason this occurs and something I should be doing differently, or is this a bug?

Thanks

-----Original Message-----
From: John Kane [mailto:jrkrideau at inbox.com] 
Sent: Thursday, July 11, 2013 10:57 AM
To: Doran, Harold; r-help at r-project.org
Subject: RE: [R] Sparse matrix no longer sparse (Matrix Package)

The message got through but not the attachment. The R help list tends to strip off attachements for security reasons.  Files of types  txt, png, & pdf should get through.

In most cases the accepted method of sending data is to use the dput() function to output a file in the console and then copy and paste the results into your email.

So for file "dat1" one would just use dput(dat1) and paste the results into an email. 

John Kane
Kingston ON Canada
-----Original Message-----
From: hdoran at air.org
Sent: Thu, 11 Jul 2013 09:53:40 +0000
To: r-help at r-project.org
Subject: [R] Sparse matrix no longer sparse (Matrix Package)

I sent this message yesterday with an attachment allowing for 
reproduction of the issue. But I think the attachment is preventing 
the message from coming through. If anyone is interested I will 
forward the attachment directly allowing for you to reproduce the 
issue I observe.

On 7/10/13 2:38 PM, "Doran, Harold" <HDoran at air.org> wrote:

I have zero'd in on what appears to be the issue. This seems to be a 
bug in Matrix, but I am not sure yet. I am attaching files that would 
allow others to replicate this with my toy data.

Notice the elements of D1 in the attached data are all integers. It 
is a sparse, diagonal matrix.

library(Matrix)
class(D1)
[1] "ddiMatrix"
attr(,"package")
[1] "Matrix"

Now, I find the inverse of the matrix A as follows:
A <- Ir + ZtZ %*% D1
A.Inv <- solve(A, Ir)
Notice now the inverse of A remains a dgCMatrix and it is relatively 
small in size, only 33424 bytes.
class(A.Inv)
[1] "dgCMatrix"
attr(,"package")
[1] "Matrix"

object.size(A.Inv)
33424 bytes

Now, if I change an element of the matrix D1 to be non-integer, D1 
still has the same class as it did before

D1[1] <- 1.2

class(D1)
[1] "ddiMatrix"
attr(,"package")
[1] "Matrix"

Now, if I use this new version of D1 in the same calculations as 
above, notice that A.Inv is no longer a dgCMatrix but instead becomes 
a dgeMatrix. It then increases from an object of size 33424 bytes to 
an object of size 2001112 bytes!

A <- Ir + ZtZ %*% D1
A.Inv <- solve(A, Ir)
class(A.Inv)
[1] "dgeMatrix"
attr(,"package")
[1] "Matrix"
object.size(A.Inv)
2001112 bytes

What I desire is that the object A.Inv remain sparse at all times and
not
become dense. But, perhaps there is a reason this change occurs that 
I don't fully understand.

I can of course coerce it back to a sparse matrix and it reduces back 
in size.
 object.size(as(A.Inv, 'sparseMatrix'))
33424 bytes

I of course recognize it requires more memory to store floating 
points than integers, but is this large increase on the order of 
magnitude that seems about right?

Is there a reason the floating point in D1 causes for A.Inv to no 
longer remain sparse?

Thank you for your help,
Harold

-----Original Message-----
From: r-help-bounces at r-project.org 
[mailto:r-help-bounces at r-project.org]
On Behalf Of Doran, Harold
Sent: Wednesday, July 10, 2013 11:42 AM
To: r-help at r-project.org
Cc: dmbates at gmail.com; 'maechler at stat.math.ethz.ch'
Subject: [R] Sparse matrix no longer sparse (Matrix Package)

I have a large function computing an iterative algorithm for fitting 
mixed linear models. Almost all code relies on functions from the 
Matrix package. I've come across an issue that I do not believe 
previously occurred in earlier versions of R or Matrix.

I have a large, sparse matrix, A as

class(A);dim(A)
[1] "dgCMatrix"
attr(,"package")
[1] "Matrix"
[1] 12312 12312

I am in a position where I must find its inverse.  I realize this is
less
than ideal, and I have two ways of doing this

A.Inv <- solve(A, Ir) or just solve(A)

Where Ir is an identity matrix with the same dimensions as A and it 
is also sparse

class(Ir)
[1] "ddiMatrix"
attr(,"package")
[1] "Matrix"

The issue, however, is that the inverse of A is converted into a 
dense matrix and this becomes a huge memory hog, causing the rest of 
the algorithm to fail. In prior versions this remained as a sparse matrix.

A.Inv[1:5, 1:5]
5 x 5 Matrix of class "dgeMatrix"
         [,1]      [,2]      [,3]      [,4]      [,5]
[1,] 0.6878713 0.0000000 0.0000000 0.0000000 0.0000000 [2,] 0.0000000
0.6718767 0.0000000 0.0000000 0.0000000 [3,] 0.0000000 0.0000000
0.5076945 0.0000000 0.0000000 [4,] 0.0000000 0.0000000 0.0000000
0.2324122 0.0000000 [5,] 0.0000000 0.0000000 0.0000000 0.0000000
0.2139975
I could coerce this matrix to become sparse such as

AA <- as(A.Inv, 'sparseMatrix')
class(AA)
[1] "dgCMatrix"
attr(,"package")
[1] "Matrix"

AA[1:5, 1:5]
5 x 5 sparse Matrix of class "dgCMatrix"

[1,] 0.6878713 .         .         .         .
[2,] .         0.6718767 .         .         .
[3,] .         .         0.5076945 .         .
[4,] .         .         .         0.2324122 .
[5,] .         .         .         .         0.2139975

But I don't think this is best.

So, my question is why is a matrix that is sparse turning into a 
dense matrix? Can I avoid that and keep it sparse without having to 
coerce it to be sparse after it is created?

Thank you very much
Harold

sessionInfo()
R version 3.0.1 (2013-05-16)
Platform: x86_64-w64-mingw32/x64 (64-bit)

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United
States.1252 [3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C 
[5] LC_TIME=English_United States.1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] lme4_0.999999-2 Matrix_1.0-12   lattice_0.20-15

loaded via a namespace (and not attached):
[1] grid_3.0.1   nlme_3.1-109 stats4_3.0.1 tools_3.0.1

[[alternative HTML version deleted]]

______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Sparse matrix no longer sparse (Matrix Package)

Thread (9 messages)