Я не знаю, является ли это оптимальным по скорости или памяти, но, вероятно, один из простых способов:
m<-apply(array(W,c(4,5000,200)),c(2,3),"%*%",r)
> m[1,1:10]==r%*%W[1:4,1:10]
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
> m[2,1:10]==r%*%W[5:8,1:10]
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE