Re: Efficient code for operating on pairwise vector diffs

[Jaroslav Hajek]

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