It is just a scaling of matrix multiplication, if you haven't recorgnize it
A=rand(300,1000);
B=rand(300,1000);
tic
C1=(A*B')/size(A,2);
toc
% or
tic
C2 = sum(pagemtimes(reshape(A,size(A,1),1,[]),reshape(B,1,size(B,1),[])),3)/size(A,2);
toc
% or
tic
C3 = sum(reshape(A,size(A,1),1,[]).*reshape(B,1,size(B,1),[]),3)/size(A,2);
toc
norm(C1-C2,'fro')
norm(C1-C3,'fro')
