1
$\begingroup$

I would like to make an Exponentially Weighted version of Correlation (similar to https://pandas.pydata.org/pandas-docs/version/0.17.0/generated/pandas.ewmcov.html).

My code is slow, how can I speed it up?


Vector Input

To do this I need to first calculate Exponentially Weighted Standard Deviations and Covariance first.

stDevEWMA[list_List, c_Real] :=  Sqrt[ExponentialMovingAverage[list^2, c]];
covarianceEWMA[list1_List, list2_List, c_Real] := ExponentialMovingAverage[list1*list2, c];

So:

correlationEWMA[list1_List, list2_List, c_Real] :=
    Last@covarianceEWMA[list1, list2, c]/( Last@stDevEWMA[list1, c]*Last@stDevEWMA[list2, c]);

Some sample data:

data = {
  {0.0025368051270555636`, 0.01974194384702299`, -0.010201971309342195`,
    0.01624296059783248`, -0.048181859416501105`, -0.02555248369410701`,
    0.011751625716882819`, 0.0009874825413342947`, 0.036892999971493534`,
    0.02918373607710656`, -0.0031430425630201153`, -0.009547744422772064`
  },
  {0.0022753767547250003`, -0.01829191934114127`, 0.023909982049822087`,
    -0.01070484309176889`, 0.04987999448021485`,  0.032412474523605406`,
    -0.015221863975058536`,0.001941816623721282`, -0.05757831949596204`,
    -0.014116594610274036`, -0.020368701580867676`, 0.0011471796777184906`
  },
  {
    0.010647727268083162`, -0.016147947259917972`, 0.011578090975155497`,
    -0.025274169235405375`, 0.05823646283591277`, 0.041215591847145294`,
    -0.02312948393856451`, -0.021664036631744432`, -0.023133482190751398`,
    -0.014041150237408373`, -0.01717626765563618`, -0.02935551962838745`
  }
};

ListLinePlot[data,PlotLegends -> Automatic]

enter image description here

The first two paths look like they are negatively correlated

correlationEWMA[data[[1]], data[[2]], 0.1]

-0.872287


When doing the calculation with two vectors as an input, I am about 20x slower than Correlation

(correlationEWMA[data[[1]], data[[2]], 0.1] // RepeatedTiming // First)/
 (Correlation[data[[1]], data[[2]]] // RepeatedTiming // First)

23.4523

Things are even worse when I adapt this to take a matrix as an input...


Matrix Input

My strategy is too loop through the upper triangular portion of the correlation matrix and then populate the remaining fields as correlation matrices are symmetric and have 1 on the diagonal

correlationMatrixEWMA[mat_List, c_Real] :=
 With[
  {
   n = Dimensions[mat][[2]]
   },
  Module[
   {
    m = ConstantArray[0, {n, n}]
    },
   Do[
    Do[
     m[[row, col]] = 
      correlationEWMA[mat[[All, row]], mat[[All, col]], c]
     ,
     {col, row + 1, n}
     ],
    {row, 1, n - 1}
    ];
   m + Transpose[m] + IdentityMatrix[n]
   ]
  ];

correlationMatrixEWMA[dataTransposed, 0.1] // MatrixForm

enter image description here

This approach is about 40x slower.

(correlationMatrixEWMA[dataTransposed, 0.1]; // RepeatedTiming // 
   First)/(Correlation[dataTransposed]; // RepeatedTiming // First)

37.3897

Is there a cleverer "Mathematica" way I can do this? Eg to make use of built in functions that are already fast (maybe WeightedData?) or do I have to go through the Compile route ?

$\endgroup$
11
  • $\begingroup$ I don't think that your definition of covarianceEWMA makes any sense because it does not center list1 and list2 in any way before ExponentialMovingAverage is applied. $\endgroup$ Jun 1, 2023 at 2:31
  • $\begingroup$ Thanks for the heads up, will check that $\endgroup$ Jun 1, 2023 at 3:29
  • $\begingroup$ covar(i) = [ a(i) * b(i) * alpha ] + [ covar(i-1) * (1-alpha) ] tallies with the definition I am using. Is this what you disagree with? I am not sure what you mean by centering the list. $\endgroup$ Jun 1, 2023 at 5:49
  • $\begingroup$ I meant that $\mathrm{cov}(X,Y)$ is $\mathrm{E} \big( (X - \mathrm{E}(X)) \, (Y - \mathrm{E}(Y)) \big)$ and not $\mathrm{E}( X, Y )$. $\endgroup$ Jun 1, 2023 at 11:41
  • $\begingroup$ I see what you mean. But are you ok with the Variance definition I use ? i.e ExponentialMovingAverage[list^2, c] $\endgroup$ Jun 1, 2023 at 13:10

1 Answer 1

1
$\begingroup$

Albeit I doubt that your implementation of correlationEWMA is correct, here is an equivalent implementation of it via Compile"

cCorrelationEWMA = Compile[{{a, _Real, 1}, {b, _Real, 1}, {q, _Real}},
   Block[{result, n, factor, p, sumaa, sumab, sumbb, ai, bi},
    n = Min[Length[a], Length[b]];
    result = Table[0., {n}];
    p = 1. - q;
    factor = q/p;
    sumaa = 0.;
    sumbb = 0.;
    sumab = 0.;
    Do[
     factor *= p;
     ai = Compile`GetElement[a, i];
     bi = Compile`GetElement[b, i];
     sumaa += ai ai factor;
     sumab += ai bi factor;
     sumbb += bi bi factor;
     , {i, n, 2, -1}];
    
    factor *= p;
    factor /= q;
    ai = Compile`GetElement[a, 1];
    bi = Compile`GetElement[b, 1];
    
    sumaa += ai ai factor;
    sumab += ai bi factor;
    sumbb += bi bi factor;
    
    sumab/Sqrt[sumaa sumbb]

    ],
   CompilationTarget -> "C",
   RuntimeOptions -> "Speed"
   ];

Timing test and correctness check:

n = 10000000;
a = RandomReal[{-1, 1}, n];
b = RandomReal[{-1, 1}, n];
q = 0.1;

result0 = correlationEWMA[a, b, q]; // RepeatedTiming // First
result1 = cCorrelationEWMA[a, b, q]; // RepeatedTiming // First
Correlation[a, b]; // RepeatedTiming // First
(result0 - result1)/result0

0.825671

0.0196983

0.0119216

-4.11029*10^-15

So it seems to provide the correct result in less that twice the amount of time spent in Correlation.

$\endgroup$
4
  • $\begingroup$ Thank you. Intuitively do you think the quickest way to build the matrix is a Compiled loop through each element? I have tried making a listable function that maps onto a vector of the upper triangular coordinates ({{1, 2}, {1, 3},...) and then building the upper triangular matrix from the results but I am not seeing much of a speed up. $\endgroup$ Jun 2, 2023 at 4:46
  • $\begingroup$ Turns out you can't Compile a ConstantArray but using your method inside the loop is very fast. +1 $\endgroup$ Jun 2, 2023 at 11:47
  • 1
    $\begingroup$ Instead of ConstantArray you can use things like Table[0., {n}], Table[0.,{m}, {n}], etc. $\endgroup$ Jun 2, 2023 at 12:25
  • $\begingroup$ Thanks, very fast now $\endgroup$ Jun 3, 2023 at 12:23

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.