I am trying for several days make this code faster. The function just takes a time series makes a Distance matrix of it and then compares components of it with Variance of that series and if component is smaller then Var, it substitutes it with 1 otherwise with 0. I would like to use it on time series of length about 10000. F.e. Dat = N[Table[Sin[x], {x, 1, 10000}]];

My function is:

Rp[data_] := (
  Clear[rr, N1, RR, DD];
  rr = Variance[data];
  N1 = Length[data];
  RR = ConstantArray[0, {N1, N1}];
  DD = DistanceMatrix[data] ^2;

  For[i = 1, i < N1, i++,
                    For[j = 1, j < N1, j++,
    If[DD[[i, j]] <= rr, RR[[i, j]] = 1; RR[[j, i]] = 1, 
     RR[[j, i]] = 0; RR[[i, j]] = 0]


I have tried several options, substitute cycle with table and compile to c, but it didn't work better. This takes on my machine cca 10 min. I can run similar code on MATLAB like in 10 sec, where is the problem please, how to speed it up ?

Thank you guys for helping me, I didn't know I can use UnitStep on whole Matrix.

But actually my whole problem is to calculate from the output Reccurence Rate-RR which is part of RQA analysis, my code is

Rqa[RP_] := (

  Lmin = 2;
  N1 = Length[RP];
  Yout = Flatten[ConstantArray[0, {1, N1}]] ;

  For[k = 2, k < N1, k++,

                  Onn = 1; While[ Onn <= N1 + 1 - k,

    If[RP[[ Onn, Onn + k - 1]] == 1, A = 1; off = 0;

                              While[ off == 0 && Onn != N1 + 1 - k, 

                               If[RP[[Onn + 1, k + Onn]] == 1,
                                          A = A + 1; Onn = Onn + 1,
                              off = 1]

                        Yout[[A]] = Yout[[A]] + 1;
    Onn = Onn + 1;
  S = 2*Yout;
  SR = 0;
  For[i = 1, i < N1, i++,  
              SR = SR + i*S[[i]];];
  RR = SR/(N1*(N1 - 1))


Which contains many loops and is also takes too much time, like 10 min for series of 10000 length = Input 10000*10000 matrix. From Michaels adviced topic I have read that I should replace For with Table ? I am afraid this wont help like in previous example. Can you please give me some advice or built in function which can help me ? I hope editing wont delete your comments.


1 Answer 1


Your problem is that you work on full matrices while the output matrix will be sparse, i.e., it contains mostly zeroes. That would be slow in any programming language. In particular, you don't need the full distance matrix. You were much better off by using Nearest.

idx = Nearest[data->Automatic,data,{∞,rr}][[All,2;;]];

should give you a list idx of lists so that idx[[i]] contains all indices j such that R[[i,j]] == 1. The matrix R is then best built with an undocumented way to create a SparseArray. The overall function Rp could look like this:

Rp[data_] := Module[{rr, N1, idx, rp, ci},
  rr = Variance[data];
  N1 = Length[data];
  idx = Nearest[data -> Automatic, data, {∞, Sqrt[rr]}][[All,2 ;;]];
  rp = Join[{0}, Accumulate[Length /@ idx]];
  ci = Partition[Join @@ idx, 1];
  SparseArray @@ {Automatic, {N1, N1}, 
    0, {1, {rp, ci}, ConstantArray[1, Length[ci]]}}

For data = N[Table[Sin[x], {x, 1, 10000}]];, Rp[data] is computed in about two and a half seconds on my machine.


Unfortunately, I observed a bit too late that the the radius to pick data is basically the standard error of the dataset. Thus, the matrix is not really sparse so that the method is only marginally faster than J.M.'s proposal.

A slightly faster way (without SparseArray and employing less mysterious syntax) is obtained with the following function:

Rp2[data_] := Module[{rr, N1, idx, R},
  rr = Variance[data];
  N1 = Length[data];
  R = ConstantArray[0, {N1, N1}];
  idx = Nearest[data -> Automatic, data, {∞, Sqrt[rr]}][[All, 2 ;;]];
  MapIndexed[{x, y} \[Function] R[[y, x]] = 1, idx];
  • $\begingroup$ Thank you very much, already Michaels solution helped me with my problem, but actually I have find out this wasnt the only thing which was slovering me, the point is I want to caluclate from previous output RR- reccurance rate, which is part of RQA analysis $\endgroup$ Mar 29, 2018 at 12:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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