Making calculation of correlation dimension faster

Question

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?

Answer 1

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

Answer 2

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