# How to improve performance of making an Exponentially Weighted Correlation Matrix?

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]


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


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 ?

• 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. Commented Jun 1, 2023 at 2:31
• Thanks for the heads up, will check that Commented Jun 1, 2023 at 3:29
• 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. Commented Jun 1, 2023 at 5:49
• 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 )$. Commented Jun 1, 2023 at 11:41
• I see what you mean. But are you ok with the Variance definition I use ? i.e ExponentialMovingAverage[list^2, c] Commented Jun 1, 2023 at 13:10

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 = CompileGetElement[a, i];
bi = CompileGetElement[b, i];
sumaa += ai ai factor;
sumab += ai bi factor;
sumbb += bi bi factor;
, {i, n, 2, -1}];

factor *= p;
factor /= q;
ai = CompileGetElement[a, 1];
bi = CompileGetElement[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.

• 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. Commented Jun 2, 2023 at 4:46
• Turns out you can't Compile a ConstantArray but using your method inside the loop is very fast. +1 Commented Jun 2, 2023 at 11:47
• Instead of ConstantArray you can use things like Table[0., {n}], Table[0.,{m}, {n}], etc. Commented Jun 2, 2023 at 12:25
• Thanks, very fast now Commented Jun 3, 2023 at 12:23