Asked by DaveF
on 18 Jan 2013

Would like to know how this code be vectorized:

X=zeros(nn); for i=1:nn for j=1:nn X(i,j)=x(i)-x(j); end end

*No products are associated with this question.*

Answer by Jan Simon
on 18 Jan 2013

Edited by Jan Simon
on 19 Jan 2013

This cannot compete with the BSXFUN approach, but it is 32% faster than the original loop:

nn = numel(x); X = zeros(nn); for i = 1:nn for j = i+1:nn a = x(i) - x(j); X(i,j) = a; X(j,i) = -a; end end

Using a temporary variable for `x(i)`, which is updated in the outer loop only, is slightly slower.

Cedric Wannaz
on 19 Jan 2013

Well then this one would be more efficient:

X = zeros(nn); for ii = 1:nn X(:,ii) = x - x(ii); end

It almost competes with bsxfun for nn=1e4.

Jan Simon
on 19 Jan 2013

@Cedric: When you post this as an answer, you get my vote.

Your version does not exploit the symmetry as the other approaches. But this modification is much slower:

X=zeros(nn); for ii = 1:nn a = x(ii:nn) - x(ii); X(ii:nn, ii) = a; X(ii, ii:nn) = -a; end

Obviously the overhead for explicit indexing is too high.

Cedric Wannaz
on 19 Jan 2013

@Jan : thank you! Knowing that I would get your vote is enough, so I'll keep that as a comment ;-)

I just profiled your last version and it is very interesting indeed, because I would never have expected such an overhead. It might even lead me to review my code in a few projects, so thank you for the comment!

Answer by Jan Simon
on 19 Jan 2013

And a C-Mex, which has the same speed as BSXFUN for a [1 x 1e3] vector under Matlab R2009a/64, Win7, MSVC 2008:

#include "mex.h"

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { double *X, *Y; mwSize i, j, n;

n = mxGetNumberOfElements(prhs[0]); X = mxGetPr(prhs[0]);

plhs[0] = mxCreateDoubleMatrix(n, n, mxREAL); Y = mxGetPr(plhs[0]);

for (i = 0; i < n; i++) { for (j = 0; j < n; j++) { *Y++ = X[i] - X[j]; } } }

Answer by Azzi Abdelmalek
on 18 Jan 2013

Edited by Azzi Abdelmalek
on 18 Jan 2013

x=[1 3 4 10 20 30]; X=cell2mat(arrayfun(@(y) y-x,x','un',0))

%or

X=bsxfun(@minus,repmat(x',1,numel(x)),x)

Cedric Wannaz
on 18 Jan 2013

Or

[X, Xt] = meshgrid(x, x) ; Xt-X

More efficient than `cell2mat/arrayfun` and less efficient than `bsxfun`. So the latter (bsxfun) should be the winner ;) But Sean's version without repmat beats them all.

Sean de Wolski
on 18 Jan 2013

Well the `for`-loops will smoke arrayfun/cell2mat like a Lambo would smoke bicycle...

`bsxfun` will be faster than the `for`-loops some of the time; meshgrid will only be faster if you need to reuse the gridded matrices.

Opportunities for recent engineering grads.

## 0 Comments