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 = 
        Length[Drop[Nearest[Drop[data, i - 1], data[[i]], {All, r}], 1]], 
        {i, 1, Length[data]}]]}, {r, 0.001, 0.011, 0.001}];

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}}, 
           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?

  • $\begingroup$ 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.) $\endgroup$ Commented Dec 26, 2017 at 0:36
  • $\begingroup$ 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! $\endgroup$
    – Michael E2
    Commented Dec 26, 2017 at 4:06
  • $\begingroup$ I believe your issue of having slower compiled code comes from the fact that Nearest isn’t compilable. $\endgroup$
    – b3m2a1
    Commented Dec 26, 2017 at 15:03
  • 1
    $\begingroup$ Related: (125543) $\endgroup$
    – corey979
    Commented Dec 26, 2017 at 20:39

2 Answers 2


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:

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

(* {1000, 3} *)

Compute the distance matrix and take the upper triangular values:

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

(* {0.051459, Null} *)

Define the correlation dimension integral function:

CDIntegral[dvals_SparseArray, n_Integer, r_?NumericQ] :=      
   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];

 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]

 (* 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"}]

enter image description here

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:

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

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


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

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


result == result2


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}];]


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


  • $\begingroup$ Nice, much faster than my implementation, especially with larger n, e.g. n=10000. $\endgroup$ Commented Dec 27, 2017 at 15:41
  • $\begingroup$ 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. $\endgroup$ Commented Dec 28, 2017 at 19:34
  • $\begingroup$ 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. $\endgroup$ Commented Dec 28, 2017 at 19:41
  • $\begingroup$ @Mark, you can use Automatic in place of "Index" in versions before 10. Re parallelization / compile, Compile might make Anton's method faster. $\endgroup$
    – kglr
    Commented Dec 28, 2017 at 20:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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