3
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ Hmm, should X=rand(500,2) instead be X=rand(n,2)? $\endgroup$ – rchilton1980 Mar 27 at 17:39
  • $\begingroup$ @rchilton1980 That's correct $\endgroup$ – EternusVia Mar 27 at 17:53
  • $\begingroup$ Are we allowed to assume that norm(X(i,:)-X(j,:),2)^2 does not overflow, I presume? $\endgroup$ – Federico Poloni Mar 28 at 23:02
  • $\begingroup$ X only has two columns, so yes. $\endgroup$ – EternusVia Mar 28 at 23:08
  • 4
    $\begingroup$ 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. $\endgroup$ – smh Mar 29 at 12:55
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Makes perfect sense. Thanks! $\endgroup$ – EternusVia Mar 27 at 19:56
1
$\begingroup$

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
$\endgroup$
  • $\begingroup$ Cool! Thank you $\endgroup$ – EternusVia Mar 28 at 23:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.