# Making calculation of correlation dimension faster

I am trying to make my calculation faster. I was hoping, when I compiled it, it would work much faster, but in the end its even slower. I don't know if I am using Compile correctly. I am open to any suggestion on how to make this computation faster.

I use the data:

DD = RandomReal[{0, 1}, {1000, 3}];


I have written the following function for calculating correlation dimension (D2):

CORRDIM[data_] := (
m =
ParallelTable[{r,
Total[
Table[
Length[Drop[Nearest[Drop[data, i - 1], data[[i]], {All, r}], 1]],
{i, 1, Length[data]}]]}, {r, 0.001, 0.011, 0.001}];
m);


CORRDIM works pretty fast, but for a large dataset it is still too slow.

CORRDIM[DD]; // AbsoluteTiming

{0.483119, Null}


I was hoping that Compile could help me speed up this calculation, especially that the options

CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True


would somehow run the code in multiple threads, but the compile version works even slower than the interpreted version.

fcc =
Compile[{{x, _Real, 2}},
Table[
{r,
Total[
Table[
Length[Drop[Nearest[Drop[x, i - 1], x[[i]], {All, r}], 1]],
{i, 1, Length[x]}]]}, {r, 0.001, 0.011, 0.001}],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True]

Map[fcc, {DD}]; // AbsoluteTiming

 {0.787973, Null}


Is it there any way to make it faster? Am I making some mistakes in the use of Compile?

• Which definition of "correlation dimension" do you use? The one I found seems to be different than what you programmed. The correlation integral in that article, $C(N,r)$, can be computed with matrix arithmetic. (Much faster than using Nearest.) Dec 26, 2017 at 0:36
• Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! Dec 26, 2017 at 4:06
• I believe your issue of having slower compiled code comes from the fact that Nearest isn’t compilable. Dec 26, 2017 at 15:03
• Related: (125543) Dec 26, 2017 at 20:39

Here is my take on this... I am following the definitions and computational steps given in the article

[1] James Theiler, "Efficient algorithm for estimating the correlation dimension from a set of discrete points", Phys. Rev. A 36, 4456 – Published 1 November 1987. DOI:https://doi.org/10.1103/PhysRevA.36.4456 .

Make data:

SeedRandom[33435]
n = 1000;
DD = Sort@RandomReal[{0, 1}, {n, 3}];
Dimensions[DD]

(* {1000, 3} *)


Compute the distance matrix and take the upper triangular values:

AbsoluteTiming[
dmat = DistanceMatrix[DD];
dvals = SparseArray[
SparseArray[UpperTriangularize[dmat, 1]]["NonzeroValues"]];
]

(* {0.051459, Null} *)


Define the correlation dimension integral function:

Clear[CDIntegral]
CDIntegral[dvals_SparseArray, n_Integer, r_?NumericQ] :=
Block[{res},
res = Total@UnitStep[r - dvals];
res*2./(n*(n - 1))
];


Note that the function CDIntegral is defined in such a way so dvals can be reused.

Compute the correlation integrals for a range of distance values.

rRange = Range[0.001, 0.011, 0.001];

AbsoluteTiming[
res =Table[CDIntegral[dvals, n, r], {r, rRange}];
]

(* {0.034141, Null} *)


This seems to be ~10 times faster, than using Nearest (on my recent MacBook Pro with \$Version 11.2.0.)

Using n=10000 the above computations are done for ~20 seconds.

Following [1] here is how the correlation dimension is estimated with the results obtained from the computations above:

lm = LinearModelFit[
Log@Select[Transpose[{rRange, res}], #[[2]] > 0 &], {1, x}, x]
lm["BestFitParameters"][[2]]

(* 2.68916 *)


Here is a plot of the found correlation integrals and the fit:

ListLogLogPlot[{Transpose[{rRange,res}], {#, Exp[lm["Function"][Log[#]]]} & /@ rRange},
Filling -> {1 -> Bottom}, PlotTheme -> "Detailed",
PlotLegends -> {"correlation integrals", "fit"}]


ClearAll[corrDIM]
corrDIM[data_] := Module[{nF = Nearest[data -> "Index"]},
Table[ {r, Tr[Table[Tr[UnitStep[nF[data[[i]], {All, r}]- i - 1]],
{i, 1, Length[data]}]]}, {r, 0.001, 0.011, 0.001}]]


Using Anton's example data:

SeedRandom[33435]
n = 1000;
DD = Sort@RandomReal[{0, 1}, {n, 3}];

First @ AbsoluteTiming[result = corrDIM[DD];]


0.044153

Comparing with OP's CORRDIM (on Wolfram Cloud where ParallelTable is not supported):

Quiet @ First @ AbsoluteTiming[result2 = CORRDIM[DD];]


0.192319

result == result2


True

Comparing with Anton's CDIntegral:

First @ AbsoluteTiming[dmat = DistanceMatrix[DD];
dvals = SparseArray[ SparseArray[UpperTriangularize[dmat, 1]]["NonzeroValues"]];
res = Table[CDIntegral[dvals, n, r], {r, rRange}];]


0.07153

2./(n*(n - 1)) result[[All, 2]] == res


True

• Nice, much faster than my implementation, especially with larger n, e.g. n=10000. Dec 27, 2017 at 15:41
• Thanks, can but I have a problem, I can not run your corrDIM function, ther is an error with "Index", what does this mean please ' Nearest[data -> "Index"] ' ? I couldnt find out. Dec 28, 2017 at 19:34
• Actually this helps me already, but originally I was trying to run computing of more of Datasets in threads, by using any of D2 estimating methods. I mean for example for all Datasets DD, DD1 DD2 ... DDn at once so it would take like the same time as for one DD. Is it possible ? I think using Compile it should work somehow, but I cant find out how. Dec 28, 2017 at 19:41
• @Mark, you can use Automatic in place of "Index" in versions before 10. Re parallelization / compile, Compile` might make Anton's method faster.
– kglr
Dec 28, 2017 at 20:07