[R] writing spdiags function for R
Moreno I. Coco
M.I.Coco at sms.ed.ac.uk
Thu Apr 12 20:10:47 CEST 2012
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.
More information about the R-help
mailing list