[R] Q. about the "solve" function in _package_ "Matrix"

[Martin Maechler]

>>>>> "yyan" == yyan liu <zhliur at yahoo.com>
>>>>>     on Wed, 22 Nov 2006 15:04:33 -0800 (PST) writes:

    yyan> Hi:
    yyan> I have some problems when I use the function "solve" function in a loop. In the following code, I have a diagonal martix "ttt" whose  elements change in every iteration in a loop. I defined a "dpoMatrix"class  before the loop so I do not need to define this class every time in the loop. The reason is to save some computing time. The code is below. The inverse of the matrix "ttt"
    yyan> should change according to the change of "ttt" in the loop. However, the values in "sigma.dpo.solve", which is the inverse of "ttt" does not change.  Can anybody figure out what is wrong with my code? 
 
Well, you should not assign to the 'x' slot in this case
since if you look at  
      str(sigma.dpo)
you see that the matrix also has "cached" its Cholesky factor,
and the direct slot assignment does not ``see the need'' for
recomputing the Cholesky factor.

Of course, one could argue this to be a bit unfortunate,
but then, you really should only directly assign to the slots of
an S4 object if you know what you are doing ... which you did
not :-) ;-)
A more extreme view would say this to be a design flaw in the "Matrix" 
*package* (not "library") that the authors have to consider.
Ideally, almost any slot reassignment of such a "Matrix" should
automatically annihilate the 'factors' slot.

  yyan>  Thanks a lot in advance!

you're welcome.
Martin Maechler, ETH Zurich

BTW, I hope that your real application does not work with
diagonal matrices, because these are much more efficiently
handled by Diagonal() and the "diagonalMatrix" class objects it
produces.

Here's your code amended  to do what you want.

## -----------------------------------------------
library(Matrix)

ttt <- diag(1,2)
str( sigma.dpo <- as(ttt, "dpoMatrix") )

## use one of these:
see.more <- TRUE
see.more <- FALSE

for(i in 1:5) {
    cat("\n-------------",i,"-----------\n")
    ttt <- diag(i,2)
    ##  assigning to 'x' slot is ``dangerous''
    ## If you really want to do this, you also need to
    ## eliminate the cashed Cholesky factor:
    sigma.dpo at x <- as.vector(ttt)
    sigma.dpo at factors <- list()
    ##
    Isigma.dpo <- solve(sigma.dpo)
    print(sigma.dpo)
    if(see.more) ## see what's going on:
        str(sigma.dpo)
    print(Isigma.dpo)
}

1
## -----------------------------------------------