This does seem like the kind of thing that there should be a faster way to compute, particularly since you are binning the results up. One approach would be to just bin the data before doing the computation, however, that loses a lot of accuracy. It does seem like there should be some moment-like approach that would allow you to bin the data before you do the computation, computing the first and second moments, or something similar and then computing the results from the binned moments. I don't know that that would work -- it's just a vague hunch. I don't know if I'll have time to try it out, but I thought I would mention it.
I
will probably now look for ways to calculate the variogram from some
random samples of my data. Thanks for the observation regarding the
square array, that would have bitten me, later.