CorrelationFunction for vectors

Question

Mathematica 9 has a new CorrelationFunction. Sadly the joy does not last long, as I can't get it to work with vectors. I would like to make a velocity correlation function. That is, given a list of n vectors I would like to calculate:

For example, this is almost instantaneous:

testdata = RandomVariate[BinormalDistribution[{10, 20}, 0.4], {10000}];
AbsoluteTiming[
x = CorrelationFunction[testdata[[All, 1]], {0, 1000}];
y = CorrelationFunction[testdata[[All, 2]], {0, 1000}];
]
ListPlot[{x, y}, PlotRange -> {0, 1}]

^{(the test data here is not the best as it is by definition uncorrelated)}

Here I would like the autocorrelation of 2D vectors evenly spaced in time. But this does not work:

r = CorrelationFunction[testdata, {0, 100}];

How can I make an autocorrelation for vectors? Or do I have to go the manual way?

ClearAll[AbsoluteAutocorrelationT, AutocorrelationT]

AbsoluteAutocorrelationT[l_, 0] := Total[(#.#) & /@ l];
AbsoluteAutocorrelationT[l_, n_] := 
  Total@Table[l[[i]].l[[i + n]], {i, Length[l] - n}];
AutocorrelationT[l_, range_] := Block[{t0, ti},
   t0 = AbsoluteAutocorrelationT[l, 0];
   ti = AbsoluteAutocorrelationT[l, #] & /@ (Range @@ range);
   (ti/t0)];

AbsoluteTiming[r = AutocorrelationT[testdata, {1, 1000}];]

Which is 200 times slower. I know I can compile this (but than I have to worry about overflows in the sum, etc.)

Answer 1

When each $v_i = (v_{i,1}, v_{i,2}, \ldots, v_{i,k})$ is a $k$-vector we may interchange the summations to obtain

$$\sum_{i=1}^{n-t} v_i\cdot v_{i+t} = \sum_{j=1}^k\sum_{i=1}^{n-t} v_{i,j} v_{i+t,j},$$

whence the desired correlation is the sum of the correlations of the components. Thus:

kernel = PDF[NormalDistribution[0, 1]] /@ Range[-5, 5, .005]; (* Induces autocorrelation *)
testdata = ListConvolve[kernel, #] & /@
  Transpose[RandomVariate[BinormalDistribution[{10, 20}, 0.4], {10^4}]];
c = Sum[ListCorrelate[v, v, {1, 1}], {v, testdata}]; // AbsoluteTiming

$\{0.0030002, \text{Null}\}$

ListLinePlot[c[[;; Length[c]/2]]/c[[1]]]

The timings scale as expected ($k=2$ shown):

Answer 2

The following is (only) six times faster:

l = RandomVariate[BinormalDistribution[{10, 20}, 0.4], {10000}];

c[l_, 0] := ListCorrelate[l, l, {1, -1}, {}, Dot, Plus, 1];
c[l_, s_List] := Module[{m = c[l, 0][[1]]}, 
   ListCorrelate[l[[;;Length@l - #]], l[[1+#;;]], {1, -1}, {}, Dot, Plus, 1] / m & /@ Range @@ s];

ListLinePlot@Table[ AbsoluteTiming[AutocorrelationT[l, {1, n}];][[1]]/
                    AbsoluteTiming[c[l, {1, n}];][[1]], {n, 30}]