Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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.

comparison graph

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"}]
share|improve this question
up vote 8 down vote accepted

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

f2[p_] := Transpose[mat[[;; -(p + 1)]]].mat[[p + 1 ;;]]
share|improve this answer
    
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

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.