Speed-up the computation of this sum of small matrices

Question

Given a matrix mat with a large number of rows (a few thousands) and a few columns (between 2 and 10), I'd like to compute the sum of the "small" matrices obtained by taking the dot product of each row (index i) and the transpose of the delayed row (index i+p):

Sum[Transpose[{mat[[i]]}].{mat[[i + p]]}, {i, 1, n - p}]

The above implementation is surely not optimal and takes quite some time to compute. How can I make it more efficient?

Example code:

n = 20000;
mat = RandomReal[{0, 5}, {n, 3}];
f[p_] := Sum[Transpose[{mat[[i]]}].{mat[[i + p]]}, {i, 1, n - p}]
Array[f, 100]; // AbsoluteTiming
(* about 15 seconds *)

Note: in signal processing,f corresponds to what is called correlation matrix function but I did not find an efficient way of using the function Correlation.

---

Edit: comparison of the answers Here is a comparison of the original solution and the two answers (log plot). It appears tom's answer is by far the most efficient.

Code used for the benchmark:

original[p_, mat_] := 
  With[{n = Length[mat]}, 
   Sum[Transpose[{mat[[i]]}].{mat[[i + p]]}, {i, 1, n - p}]];
drbelisarius[p_, mat_] := 
  With[{n = Length[mat]}, 
   Inner[Times, mat[[1 ;; n - p]], mat[[p + 1 ;; n]], Plus, 1]];
tom[p_, mat_] := Transpose[mat[[;; -(p + 1)]]].mat[[p + 1 ;;]];
compareStep[p_, mat_] := 
 Map[AbsoluteTiming[#[p, mat];] &, {original, drbelisarius, 
    tom}][[All, 1]]
SeedRandom[6666]
compare[n_] := Block[{mat = RandomReal[{0, 5}, {n, 3}]},
  Array[compareStep[#, mat] &, n - 1] // Total]
tab = Array[compare, 300][[2 ;; -1]];
ListLogPlot[Transpose@tab, AxesLabel -> {"n", "time"}, 
 PlotLegends -> {"original", "drbelisarius", "tom"}]

Answer 1

You can write the summation and the outer product as a matrix multiplication:

f2[p_] := Transpose[mat[[;; -(p + 1)]]].mat[[p + 1 ;;]]

Answer 2

ff[p_] := Inner[Times, mat[[1 ;; n - p]], mat[[p + 1 ;; n]], Plus, 1]

Seems 3x faster