[R] diag()<- in Matrix?

[Martin Maechler]

>>>>> "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