[R] writing spdiags function for R

[Moreno I. Coco]

Dear R-list,

I am in the process of translating a long function written in Matlab
into R (mainly because I am a big of fan of R, and folks will not
have to pay to use it :). In the translation of this function
I got stack because they use spdiags, which, as far as I can tell
it is not available in R. I have explored the Matrix package, from
which I borrowed some of the functions (e.g., sparseMatrix), but
I could not actually find an equivalent to spdiags (please, let
me know if it is there somewhere).

So, I have written my own spdiags function (below); following
also a suggestion in an old, and perhaps unique post, about
this issue.

It works only for square matrices (that's my need), however I
have a couple of issues, mainly related to computational
efficiency:

1) if I use it with a sparseMatrix, it throws a tedious warning
"<sparse>[ <logic> ] : .M.sub.i.logical() maybe inefficient";
can I suppress this warning somehow, this is slowing the computation
very radically;

2) I can go around this problem by translating a sparseMatrix back
into a logical matrix before I run spdiags on it. However, the loop
gets very slow for large matrices (e.g., 2000x2000), which is the
kind of matrices I have to handle. If you look in the code,
I have placed a system.time() where the code is slowing down, and
it takes about:

user  system elapsed
    0.28    0.05    0.33

to complete an iteration...thus, I was wondering whether there is
a more efficient way to do what I am doing...also, if you spot
other places where the function could be optimized I would be
very glad to hear it!

thank you very much in advance for your kind help,

Best,

Moreno

#######################################################################

## it works only for square matrices
## it could work with sparse matrices but it spits a tedious warning
## it is definitely inefficient compared to the original matlab code

## choose below different matrices to test the function.
# r = c(2,3,5,5); c = c(2,1,4,5)
# A = sparseMatrix(r, c)
# A = replicate(1000, rnorm(1000) )
# A = rbind(c(1,2,3),c(2,3,4),c(3,4,5))

spdiags = function(A){

     # Find all nonzero diagonals
     i = rep(seq(1, nrow(A),1),nrow(A));
     j = sort(i);
     d = sort(j-i);

       # d = d(find(diff([-inf; d]))); ## from Matlab ...
       # d = c(d[which(diff(d) == 1)], d[length(d)] ) ## this emulate 
above but needs to stick in last element

     d = unique(d); ##this should work just fine and it is simpler
     p = length(d); ##the nr. col of the new matrix
     m = nrow(A); n = ncol(A);

     B = matrix(0, nrow = min(c(m,n)), ncol = p);

   for (k in 1:p){
      # print(k)
      cl = vector();

       if (m >= n){
          i = max(1, 1+d[k]):min(n, m+d[k]);
       } else { i = max(1, 1-d[k]):min(m,n-d[k]); }

       system.time(
       if (length(i) != 0){
          B[i,k] = A[ col(A) == row (A) - d[k]]
       } )
}

return (list( B = B, d = d) )

}




-- 
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.