[R] diag()<- in Matrix?
Martin Maechler
maechler at stat.math.ethz.ch
Mon Nov 6 15:58:38 CET 2006
>>>>> "Jerome" == Jerome Goudet <jerome.goudet at unil.ch>
>>>>> on Sun, 05 Nov 2006 16:08:33 +0100 writes:
Jerome> Dear all,
Jerome> I am trying to use the Matrix package to do some calculations on rather
Jerome> large and sparse matrices. An example of such a matrix is given below.
Jerome> mig<-0.2
Jerome> side<-10
Jerome> np<-side^2
Jerome> mig.mat<-matrix(0,np,np)
Jerome> diag(mig.mat[1:(np-side),(side+1):np])<-mig/4
Jerome> diag(mig.mat[(side+1):np,1:(np-side)])<-mig/4
Jerome> diag(mig.mat[-np,-1])<-mig/4
Jerome> nomig<-which(1:(side^2-1)%%side==0)
Jerome> diag(mig.mat[-np,-1])[nomig]<-0
Jerome> diag(mig.mat[-1,-np])<-mig/4
Jerome> diag(mig.mat[-1,-np])[nomig]<-0
Jerome> diag(mig.mat)<-1-apply(mig.mat,2,sum)
Jerome> Here, side is set to 10, but ideally, I'd like it to be 500.
Jerome> I tried to simply replace matrix by Matrix in the code above, but
Jerome> obviously, it did not work...
^^^^^^^^^
Jerome> The error came from the line below:
Jerome> diag(mig.mat[1:(np-side),(side+1):np])<-mig/4
Jerome> the diag part is ok, but it seems that I cannot assign values to this
Jerome> vector.
Yes, "diag<-", the function which is called for this, does not
work for 'Matrix' matrices, actually for two reasons :
1) it has check which is not appropriate for 'Matrix' matrices
2) it is based on
M[ cbind(i,j) ] <- value
the famous not-well-enough-known other way of indexing in S
and R,
but this way of indexing and index-assignment has been on our
TODO list, but has not yet been implemented for 'Matrix'
matrices.
Your use-case is a good incentive for
Doug Bates and me (the 'Matrix' authors) to implement it -- for
a next version of the Matrix package
{ BTW: your version 0.9975-3 is outdated anyway,
and you *should* update.packages()
}
Because of the above, I'll probably also be able to help you in
a short while with a workaround before this works
out-of-the-box,
but please do upgrade to the latest released version of
Matrix (0.9975-5), in any case.
Jerome> I would be very glad for any help on this, or on pointers on how to
Jerome> construct Matrices.
Jerome> As I am not on the list, please could you answer to
Jerome> jerome.goudet at unil.ch?
Jerome> Many thanks in advance.
you're welcome:
Martin Maechler, ETH Zurich
Jerome> R 2.4.0
Jerome> Windows XP
Jerome> Matrix package version 0.9975-3
More information about the R-help
mailing list