[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Efficient code for operating on pairwise vector diffs
From: |
Jaroslav Hajek |
Subject: |
Re: Efficient code for operating on pairwise vector diffs |
Date: |
Mon, 2 Nov 2009 15:00:08 +0100 |
On Mon, Nov 2, 2009 at 2:21 PM, Joseph Wakeling
<address@hidden> wrote:
> Hello all,
>
> Suppose that I have a set of N d-dimensional vectors, which I guess
> could be represented as either a matrix,
>
> X = [X1 X2 ... XN]
>
> or a cell:
>
> X = {X1 X2 ... XN}
>
> I'm trying to implement a function which consists of the following steps:
>
> (i) calculate the pairwise vector differences Dij = Xi - Xj
>
> (ii) perform some function on each of the vector differences,
> Fij = F(Dij)
>
> (iii) Construct the averages aveFi = mean of Fij calculated over
> the subscript j.
>
> My guess is that the best way to do this is to represent the vector
> differences in a NxN cell whose (i,j)'th element is the d-dimensional
> vector X(:,i)-X(:,j). Then I can call cellfun to apply the function F
> and finally the averages can be achieved by some use of mean(cell2mat())
> or whatever.
>
> Now, the challenge is how to generate the vector differences
> efficiently. I can do,
>
> for i=1:N
> for j=1:N
> D{i,j} = X(:,i) - X(:,j)
> end
> end
>
> ... but that seems very inefficient, which is exactly what I'm trying to
> avoid(*). Can anyone advise an efficient way to carry out the above?
>
> Thanks & best wishes,
>
> -- Joe
If X is a DxN matrix, you can get a dxNxN array of pairwise
differences using, for instance using spread indexing
D = X(:,:,ones (1, N)) - reshape (X, d, 1, N)(:,ones (1, N), :);
or using bsxfun
D = bsxfun (@minus, X, reshape (X, d, 1, N));
the relative efficiency will depend on your version of Octave (bsxfun
wins if you use the bleeding edge sources). Then,
D = squeeze (nu2mcell (D, 1));
converts this to an NxN cell array of dx1 vectors. Note that unless
your mapper function is fast enough (built-in?), this isn't going to
win you much, because your function will still be called NxN times.
If F is a fixed mapping and relatively simple, there may be a better
way to do the trick in a few calls. See also
http://octave.sourceforge.net/doc/f/pdist.html
> (*) Motivation for this whole query: I'm editing some MATLAB code by a
> colleague which we're going to release some time soon(**). His code
> works but is not very efficiently written from a MATLAB/Octave point of
> view, and so he's used C extensions for some parts.
>
> I think that the code can be made efficient without needing to use C,
> which would be very nice as it would allow it to be run in Octave as
> well as MATLAB. But my implementation is foundering on implementing
> the above code without costly for loops ...
>
See the above comment. Note that Octave also has a C++ extension as an option.
> (**) In answer to 'When?' and 'Why can't we see the actual code now?':
> it comes down to academic confidentiality. If it was just my work I
> would happily share, but I can't presume on the part of my co-workers.
Code is always better if you want good answers. You can only extract a
small (but working!) skeleton. Unless your co-workers are idiots, they
won't bark at you for showing a few lines, especially if you're trying
to improve your common work.
hth
--
RNDr. Jaroslav Hajek
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz