# Efficient covariance matrix calculation MATLAB (every combination of rows from data)

My friend in the statistics department asked me how to do the following calculation efficiently.

Suppose we have data $$X\in\mathbb{R}^{N\times 2}$$. He needs to do the following calculation: $$C_{i,j}=\sigma^2\exp\left(-\frac{||X_i-X_j||_2}{\theta}\right)$$ where $$X_i$$ is a tuple with two elements and $$\sigma,\theta$$ are scalars. So $$C$$ is a $$N\times N$$ matrix computed by taking every possible combination of rows from $$X$$.

I fancy myself a MATLAB aficionado, but I could not figure out how to efficiently compute this using the built in MATLAB matrix operations and functions. The fastest implementation I have uses the standard for-loops:

clear all; close all; clc;

n=10;
X=rand(n,2);
sigma=1;
theta=1;

tic
C=zeros(n,n);
for i=1:n
for j=1:n
C(i,j)=sigma*sigma*exp(-norm(X(i,:)-X(j,:),2)/theta);
end
end
toc


Is there a faster way?

• Hmm, should X=rand(500,2) instead be X=rand(n,2)? – rchilton1980 Mar 27 at 17:39
• @rchilton1980 That's correct – EternusVia Mar 27 at 17:53
• Are we allowed to assume that norm(X(i,:)-X(j,:),2)^2 does not overflow, I presume? – Federico Poloni Mar 28 at 23:02
• X only has two columns, so yes. – EternusVia Mar 28 at 23:08
• The Gaussian kernel matrix is known to have low-rank structure. For large $N$, e.g. $N > 5000$ or so, it may be useful to use low-rank approximation methods such as those used in H-matrices (example matlab code: github.com/marianona/Hmatrix) to calculate a compressed representation in $\mathcal{O}(N)$ time and requiring only $\mathcal{O}(N)$ memory. Such a representation is very efficient and can be arbitrarily accurate for this particular matrix. – smh Mar 29 at 12:55

Try the snippet below:

clear all;
close all;

n=500;
rand('seed',0)
X=rand(n,2);
sigma=1;
theta=1;

% Original method, A
tic
Ca=zeros(n,n);
% DX1a = zeros(n,n);
% DX2a = zeros(n,n);
for i=1:n
for j=1:n
Ca(i,j)=sigma*sigma*exp(-norm(X(i,:)-X(j,:),2)/theta);
% DX1a(i,j) = X(i,1) - X(j,1);
% DX2a(i,j) = X(i,2) - X(j,2);
end
end
toc

% Alternative method, B
tic
X1 = X(:,1);
X2 = X(:,2);
[MX1, MX2] = meshgrid(X1,X2);
DX1b = MX1'-MX1;
DX2b = MX2-MX2';
Cb = sigma*sigma*exp(-sqrt(DX1b.^2 + DX2b.^2)/theta);
toc

C_error = norm(Ca-Cb,'fro')

% DX1_error = norm(DX1b-DX1a,'fro')
% DX2_error = norm(DX2b-DX2a,'fro')


The trick here is using meshgrid to spill/splay the coordinates into a pair of NxN arrays (here MX1 for coordinate 1, and MX2 for coordinate 2) so that you can compute the coordinate-difference matrices, DX1 and DX2, all in one go. Then it's just a few elementwise ops (squares, sqrts, exp's) to form C. Although the original method didn't compute DX1 or DX2, I put some commented-out code in there, that you can turn back on if you wish to probe more deeply. Method B is considerably faster, on my box anyway:

Elapsed time is 7.594 seconds.
Elapsed time is 0.0380321 seconds.
C_error =   1.9128e-014


I encourage you to test some more. The meshgrid() function is a little idiosyncratic, so double check on your datasets / larger program just to make sure.

• Makes perfect sense. Thanks! – EternusVia Mar 27 at 19:56

A slightly faster variant of @rchilton1980's method that uses singleton expansion rather than meshgrid:

Cc = sigma * sigma * exp(-sqrt((X(:,1) - X(:,1)').^2 + (X(:,2) - X(:,2)').^2) / theta);


And yet another variant that uses vecnorm (introduced in R2017b). Its speed is on par with method B, but I guess it's going to be the most efficient solution when $$X$$ is $$n\times k$$ with $$k>2$$.

Cd = sigma * sigma * exp(-vecnorm(reshape(X, [500,1,2]) - reshape(X, [1,500,2]), 2, 3)/theta);


On my machine:

Elapsed time is 0.575712 seconds. % Method A
Elapsed time is 0.013994 seconds. % Method B
Elapsed time is 0.006800 seconds. % Method C
Elapsed time is 0.016893 seconds. % Method D
Cb_error =
1.3146e-14
Cc_error =
1.3146e-14
Cd_error =
0

• Cool! Thank you – EternusVia Mar 28 at 23:09