1
$\begingroup$

I have a simple module that I want to compile to potentially speed up.

normCrossCorr2[left0_, right0_] :=
  Module[{left = left0, right = right0, dim, leftbar, rightbar, upper,
     lower, xcorr},
   
   left = Flatten[left];
   right = Flatten[right];
   dim = Length[right];
   leftbar = Mean[left];
   rightbar = Mean[right];
   
   upper = Total[(left - leftbar)*(right - rightbar)];
   lower = Sqrt[Total[(left - leftbar)^2]*Total[(right - rightbar)^2]];
   xcorr = upper/lower
   ];

I can compile it - it appears to compile fine.

cf = Compile[{left0, right0}, normCrossCorr2[left0_, right0_]]

But when testing it with the following code, the original module runs just fine, but the compiled module errors (see image below).

a = RandomReal[{0, 1}, {10, 10}];
b = RandomReal[{0, 1}, {10, 10}];
normCrossCorr2[a, b]
cf[a, b]

This is the first time I have used Compile, so I am sure I am just being dumb.

The error is related to machine size real numbers, but all of the input is machine size and real. I therefore assume one of the inner variables created is not?, but I am not sure why.

I have had a good look at the documentation and online, but I have still not been able to solve this issue.

Any help would be appreciated.

Also any pointer to speeding up my Module would be great. The example here is toy, the real example much more computational intensive as a and b are much bigger. The Module will also hopefully run in a ParallelTable.

I have tried to be as comprehensive as possible, without providing needless detail. Please let me know if anything is unclear.

Note: the image of the output is much bigger, but it is just more of the same as far as I can see, so I have cut it short for purposes of legibility.

Thank you.

Image of error

$\endgroup$
4
  • $\begingroup$ You have to remove the pattern bars (_) in Compile[{left0, right0}, normCrossCorr2[left0_, right0_]] $\endgroup$ Apr 7, 2023 at 13:22
  • $\begingroup$ Moreover, in the definition normCrossCorr2 left0 and right0 seem to be some multidimensional arrays. Of course, you have to Compile tell about this by specifying the types of its arguments. See the documention of Compile. $\endgroup$ Apr 7, 2023 at 13:25
  • $\begingroup$ Thanks. Removing the bars helps with some of the errors, but not the machine size real error. The input arguments are each grids of numbers i.e. rows and cols. I think I am missing how to specify this. $\endgroup$
    – flyingmind
    Apr 7, 2023 at 13:34
  • 1
    $\begingroup$ You could simply use Correlation[Flatten[left0], Flatten[right0]]; Correlation is already compiled and should be quite efficient. $\endgroup$ Apr 7, 2023 at 13:38

1 Answer 1

2
$\begingroup$

This is the most economical way I can imagine to implement this in Mathematica without recoursing to LibraryLink. The idea is to have a single double loop that computes moments up to degree 2 and to compute the correlation from that. However, beware of the fact that this may lead to quite some rounding errors if the inputs X and Y are very far from being centered.

cf = Compile[{{X, _Real, 2}, {Y, _Real, 2}},
   Module[{m, n, div, Xij, Yij, XMean, YMean, CovXY, CovXX, CovYY},
    {m, n} = Dimensions[X];
    XMean = 0.;
    YMean = 0.;
    CovXY = 0.;
    CovXX = 0.;
    CovYY = 0.;
    Do[
     Xij = Compile`GetElement[X, i, j];
     Yij = Compile`GetElement[Y, i, j];
     XMean += Xij;
     YMean += Yij;
     CovXX += Xij Xij;
     CovXY += Xij Yij;
     CovYY += Yij Yij;
     , {i, 1, m}, {j, 1, n}];
    
    div = 1./N[m n];
    XMean *= div;
    YMean *= div;
    CovXY = CovXY div - XMean YMean;
    CovXX = CovXX div - XMean XMean;
    CovYY = CovYY div - YMean YMean;
    
    CovXY/Sqrt[CovXX CovYY]
    ],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];

This is already parallelized when threaded over lists of matrices. Its speed advantage over ParallelTable and Correlation depends heavily on the size of the inputs. Typically, ParallelTable has quite a high overhead so that it often pays off only for large enough jobs. Anyways, here a brief test example:

k = 80;
m = 1000;
n = 1000;
X = RandomReal[{-1, 1}, {k, m, n}];
Y = RandomReal[{-1, 1}, {k, m, n}];

r1 = ParallelTable[
     Correlation[Flatten[X[[i]]], Flatten[Y[[i]]]],
     {i, 1, Length[a]},
     Method -> "CoarsestGrained"
     ]; // RepeatedTiming // First
r2 = cf[X, Y]; // RepeatedTiming // First

Max[Abs[r1 - r2]]

0.155226

0.0306176

2.09468*10^-16

A version that is a bit less susceptable to rounding errors is this one. It takes 33% longer, but that is probably worth it.

cfRobust = Compile[{{X, _Real, 2}, {Y, _Real, 2}},
   Module[{m, n, div, Xij, Yij, XMean, YMean, CovXY, CovXX, CovYY},
    {m, n} = Dimensions[X];
    XMean = 0.;
    YMean = 0.;
    
    Do[
     Xij = Compile`GetElement[X, i, j];
     Yij = Compile`GetElement[Y, i, j];
     XMean += Xij;
     YMean += Yij;
     , {i, 1, m}, {j, 1, n}];
    
    div = 1./N[m n];
    XMean *= div;
    YMean *= div;
    
    CovXY = 0.;
    CovXX = 0.;
    CovYY = 0.;
    
    Do[
     Xij = Compile`GetElement[X, i, j] - XMean;
     Yij = Compile`GetElement[Y, i, j] - YMean;
     CovXX += Xij Xij;
     CovXY += Xij Yij;
     CovYY += Yij Yij;
     , {i, 1, m}, {j, 1, n}];
    
    CovXY = CovXY div;
    CovXX = CovXX div;
    CovYY = CovYY div;
    
    CovXY/Sqrt[CovXX CovYY]
    ],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];

Note also that Total, Mean, Covariance, Correlation, etc. have probably implemented clever means to protect against catastrophic cancellation. So if you have reason to believe to get into trouble, better use those ones. But then there is little that can be obtained from putting them into a Compile as each of these functions quite certainly has already a compiled backend.

$\endgroup$
4
  • $\begingroup$ Wonderful! Thank-you, I did to expect so much. I will have a play with this tonight/tomorrow. $\endgroup$
    – flyingmind
    Apr 7, 2023 at 15:56
  • $\begingroup$ You're welcome. $\endgroup$ Apr 8, 2023 at 7:41
  • $\begingroup$ A quick question. I notice that this is coding in what I would call a "C-like" way (and compiled to C) i.e. loops. Is this general best advice? For example, in "plain" Mathematica my understanding is that loops are not necessarily good (normally hidden via Table or Map et al). Basically: when aiming to compile code should one program in Mathematica in a "Mathematica-like" way, or program in Mathematica in a "C-like" way (or the language you are aiming to compile for)? $\endgroup$
    – flyingmind
    Apr 8, 2023 at 11:53
  • 1
    $\begingroup$ Within Compile do it C-like. The point of Compile is to use the advantages of C (strong typing, fast loops, compile time optimizations, etc). The thinner the abstraction layer, the more likely the C-code generate will generate code that the compiler understands and is able to optimize $\endgroup$ Apr 8, 2023 at 22:24

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.