# Can this Sample Entropy (SampEn) calculation be improved?

I am interested in calculating various entropy measures of time series, and was surprised to find that none of them are included with Mathematica. Sample Entropy (SampEn) is a fairly common technique for this purpose, so I've written a module to perform the calculation, but I'm surprised at how slow it is (especially compared to my Python implementation).

SampEn[TheData_,TheWindowSize_,TheThreshold_]:=
Module[{Num=Length[TheData],NumWindows,TheValueArray1,TheValueArray2,TheDifferences1,,TheDifferences2},
NumWindows=Num-TheWindowSize+1;
TheValueArray1=Table[Table[TheData[[j]],{j,i,i+(TheWindowSize-1)}],{i,NumWindows}];
TheValueArray2=Table[Table[TheData[[j]],{j,i,i+(TheWindowSize)}],{i,NumWindows-1}];
TheDifferences1=Total[Table[Total[Table[If[ChessboardDistance[i,j]<TheThreshold,1,0],{j,TheValueArray1}]]-1,{i,TheValueArray1}]];
TheDifferences2=Total[Table[Total[Table[If[ChessboardDistance[i,j]<TheThreshold,1,0],{j,TheValueArray2}]]-1,{i,TheValueArray2}]];
-Log[N@(TheDifferences2/TheDifferences1)]];
RandomWalkingData=Accumulate[Table[1.0*RandomVariate[NormalDistribution[0,.4]],{t,1,1000,1}]];
m=2;
r=0.2*StandardDeviation[RandomWalkingData];
Timing[SampEn[RandomWalkingData,m,r]]


This takes 40.04s for 1000 datapoints, but with 5000 points it already takes 989.608s. I also tried a variant that combined the Table and Total operations, but it wasn't any faster.

TheDifferences1=Total[Table[If[ChessboardDistance[i,j]<TheThreshold,1,0],{j,TheValueArray1},{i,TheValueArray1}],2]-Length[TheValueArray1];


Are there any shortcuts or built-in functions that can be used to speed up this calculation? How can this be made more efficient?

kglr did already the major part of the work by boiling down things to Nearest. I have only two further remarks:

1. NearestFunctions thread over lists, enabling parallelization.

2. The bandwidth limited (and thus expensive) list operations like Rest and Flatten can be avoided when it is just about counting the number of elements.

The resulting code looks like this:

ClearAll[sampleEnt2]
sampleEnt2[data_, ws_, t_] :=  Module[{nF1, nF2, diff1, diff2, va1, va2},
va1 = Partition[data, ws, 1];
va2 = Partition[data, 1 + ws, 1];
nF1 = Nearest[va1 -> Automatic, DistanceFunction -> ChessboardDistance];
nF2 = Nearest[va2 -> Automatic, DistanceFunction -> ChessboardDistance];
diff1 = Total[(Length /@ nF1[va1, {All, t}])] - Length[va1];
diff2 = Total[(Length /@ nF2[va2, {All, t}])] - Length[va2];
-Log[N@(diff2/diff1)]
]


This is how it performs on my quad core machine:

SeedRandom
RandomWalkingData =
Accumulate[1.0 RandomVariate[NormalDistribution[0, .4], 1000]];
m = 2;
r = 0.2*StandardDeviation[RandomWalkingData];
SampEn[RandomWalkingData, m, r] // AbsoluteTiming
sampleEnt[RandomWalkingData, m, r] // AbsoluteTiming
sampleEnt2[RandomWalkingData, m, r] // AbsoluteTiming


{36.0335, 0.120787}

{0.02856, 0.120787}

{0.006953, 0.120787}

• +1, great points. – kglr Dec 27 '17 at 7:39
ClearAll[sampleEnt]
sampleEnt[data_, ws_, t_] :=   Module[{nF1, nF2,  diff1,  diff2,
va1 = Partition[data, ws, 1], va2 = Partition[data, 1 + ws, 1]},
nF1 = Nearest[va1 ->"Index", DistanceFunction -> ChessboardDistance];
nF2 = Nearest[va2 -> "Index", DistanceFunction -> ChessboardDistance];
diff1 = Length[Flatten[Rest[nF1[#, {All, t}]]& /@ va1]];
diff2 = Length[Flatten[Rest[nF2[#, {All, t}]]& /@ va2]];
-Log[N @ (diff2 / diff1)]]


Example:

SeedRandom
ss = 500;
RandomWalkingData = Accumulate[Table[1.0 RandomVariate[NormalDistribution[0, .4]],
{t, 1, ss, 1}]];
m = 2;
r = 0.2*StandardDeviation[RandomWalkingData];
Timing[sampleEnt[RandomWalkingData, m, r]]


{0.007999, 0.482223}

This is about 1500X speed-up versus OP's SampEn on the same input:

Timing[SampEn[RandomWalkingData,m,r]]


{11.5612, 0.482223}

For a larger input size (ss = 1000) we get

Timing[sampleEnt[RandomWalkingData, m, r]]


{0.029995, 0.26345439758666905}

Timing[SampEn[RandomWalkingData, m, r]]


{46.099992, 0.26345439758666905}

• Looks good, but I get an error {.....{0.9616,0.9613,0.9716},{0.9613,0.9716,0.962},<<4299>>}->Index is neither a list of real points nor a valid list of rules. using v10.4, and actually I don't understand what ->"Index"' means here. – Aaron Bramson Dec 27 '17 at 7:30
• Aaron, try Nearest[va1 -> Automatic, DistanceFunction -> ChessboardDistance], instead. The syntax with Index was introduced around version 11.2. Both variants mean that only the positions should be looked up by Nearest`, not the points themselves. – Henrik Schumacher Dec 27 '17 at 7:31