Re: Difference of columns

[Jordi Gutiérrez Hermoso]

On Thu, 2013-08-29 at 13:30 +0200, Ismael Núñez-Riboni wrote:
> >> Suppose that we have N columns in a matrix, how we can perform a
> >> program that take into account all possible differences between these
> >> N columns.
> 
> You need a double loop, I see no other (easy) way, if your matrix is A:

I wish broadcasting were more well-known, because it is supposed to be
precisely this easy way. It is precisely why I enabled it since Octave
3.6, so that people would use it more:

    reshape (A - permute (A, [1,3,2]),
            rows (A), columns (A)^2)

This gives you both colI - colJ and colJ - colI, and also gives you a
bunch of colK - colK zero columns.

To remove them:

    reshape (A - permute (A, [1,3,2]),
             rows (A), columns (A)^2)(:, find (triu (ones(columns(A)),
                                                     1)));

I know what you're going to think. This is way too complicated, we
should be writing double loops instead. But my code above is quite
obvious after you're used to broadcasting. Suppose for now that all of
your columns were of length one and you wanted to get all pairwise
differences between your columns, i.e. between your scalars. This is
just a subtraction table:

    v = [1, 2, 3, 4, 5];
    v - v' ## This is a subtraction table.

You can still do a subtraction table if you have columns instead of
scalars. In this case, one way to think about it is that each row in
your columns becomes a stacked subtraction table. The stacked
subtraction table is thus a 3d array. You use permute to orient the
matrix in the right way 3-dimensionally and you subtract as usual,
getting a stacked 3d subtraction table:

    table = A - permute (A, [1,3,2]);

Then you flatten this table back into a matrix:

    table = reshape (A, rows (A), columns (A)^2);

However, you really only want the top triangle of this subtraction
table, so that you omit the zero diagonal and the lower triangle which
is just the negative of the top triangle. For this, you do an indexing
trick with triu and find to get the indices of this upper triangle:

    idx = triu (ones ( columns (A)), 1);
    idx = find (A); ## Flatten the index, to match flattened table

Finally, we use this index to get the columns of the flattened
subtraction table that we want:

    table = table(:, idx);

What I wrote above is the same but without intermediate variables.

HTH,
- Jordi G. H.