Speed-up the computation of this sum of small matrices

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


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"}]

-

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

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

-
Your code is much more efficient (see edit in OP). I will accept it if no (even) better solution come. – anderstood Mar 3 at 14:50
@anderstood, dot products are quite efficient; I doubt that this solution here can be beaten. – J. M. Mar 3 at 15:17
ff[p_] := Inner[Times, mat[[1 ;; n - p]], mat[[p + 1 ;; n]], Plus, 1]


Seems 3x faster

-