A better way to do antidiagonal matrix vectorization?

Say I want to reduce a MxN matrix to a column vector, consisting of all antidiagonal vectors. Consider the 5x5 matrix:

A = reshape(1:25,[5 5]) % input

A = 5×5

1 6 11 16 21 2 7 12 17 22 3 8 13 18 23 4 9 14 19 24 5 10 15 20 25

B = [1 6 2 11 7 3 16 12 8 4 21 17 13 9 5 22 18 14 10 23 19 15 24 20 25].' % intended output

B = 25×1

1 6 2 11 7 3 16 12 8 4

In actuality, I want this to also work on 3D arrays by simply doing the same operation on all pages, resulting in a (M*N)x(numpages) matrix. I should point out that all these examples are square, but the solution should work generally for any rectangular page geometry.

I came up with four solutions of varying clumsiness and inelegance. The first two approaches use diag(); the latter two use circshift().

% build 1000x1000x3 test array
N = 1000;
s = [N N];
A0 = reshape(1:prod(s),s);
A0 = repmat(A0,[1 1 3]);

% just call diag() a bazillion times for each channel
timeit(@() advectorize1(A0))

ans = 0.0423

% use an index array and only call diag() one bazillion times
timeit(@() advectorize2(A0))

ans = 0.0519

% build index array, nan pad, vector circshift
timeit(@() advectorize3(A0))

ans = 0.0831

% can work on the image itself, but generally the slowest
% no index array, but pushes around 3x as much data during circulation
timeit(@() advectorize4(A0))

ans = 0.1761

function B = advectorize1(A)
    A = flipud(permute(A,[2 1 3]));
    s = size(A);
    B = zeros(prod(s(1:2)),s(3));
    for c = 1:s(3)
        Ac = A(:,:,c);
        ipr = 0;
        for k = -(s(1)-1):(s(2)-1)
            d = diag(Ac,k);
            inx = ipr + numel(d);
            B(ipr+1:inx,c) = d;
            ipr = inx;
        end
    end
end

function B = advectorize2(A)
    A = flipud(permute(A,[2 1 3]));
    s = size(A);
    idx0 = reshape(1:prod(s(1:2)),s(1:2));
    idx = zeros(prod(s(1:2)),1);
    ipr = 0;
    for k = -(s(1)-1):(s(2)-1)
        d = diag(idx0,k);
        inx = ipr + numel(d);
        idx(ipr+1:inx) = d;
        ipr = inx;
    end
    B = A((idx + prod(s(1:2))*(0:s(3)-1)));
    B = reshape(B,[],s(3));
end

function B = advectorize3(A)
    s = size(A);
    idx = reshape(1:prod(s(1:2)),s(1:2));
    idx = [idx nan(s(1),s(1)-1)];
    for m = 2:s(1)
	    idx(m,:) = circshift(idx(m,:),m-1,2);
    end
    idx = idx(~isnan(idx)); % mask & vectorize
    B = A((idx + prod(s(1:2))*(0:s(3)-1)));
    B = reshape(B,[],s(3));
end

function B = advectorize4(A)
    A = double(A); % isn't integer-compatible
    s = size(A);
    A = [A nan(s(1),s(1)-1,s(3))];
    for m = 2:s(1)
	    A(m,:,:) = circshift(A(m,:,:),m-1,2);
    end
    B = A(~isnan(A));
    B = reshape(B,[],s(3));
end

The four methods are listed in order of decreasing speed for large (>> 200x200x3) 3D arrays. For small arrays, the order is reversed, though the shorter execution times make that case relatively unimportant.

So far, method 1 is the best. It isn't problematically slow, but It seems like there would surely be something more elegant than using diag() a whole bunch of times -- maybe even something both elegant and faster? Am I asking too much?

Any ideas?

Answer 1

[Paul]

The approach in test1() isn't faster, but perhaps it meets the criterion for more elegance. My tests indicate that most of the time is spent on the sort. So if your workflow allows you to hoist the computation of ind out of the function, e.g., if calling this function in a loop with different data but all of the same dimenstion, then test1() reduces to fairly simple operations that I think would be quite fast.

% verify operation for some simple cases
A = [1 3 5;40 50 60;70 80 90]  % square

A = 3×3

1 3 5 40 50 60 70 80 90

test1(A)

ans = 9×1

1 3 40 5 50 70 60 80 90

A = [1 3 5;40 50 60]  % wide

A = 2×3

1 3 5 40 50 60

test1(A)

ans = 6×1

1 3 40 5 50 60

A = [1 3 5;40 50 60].'  % tall

A = 3×2

1 40 3 50 5 60

test1(A)

ans = 6×1

1 40 3 50 5 60

N = 1000;
s = [N N];
A0 = reshape(1:prod(s),s);
A0 = repmat(A0,[1 1 3]);

B1 = advectorize1(A0);
B2 = test1(A0);
isequal(B1,B2)  % check

ans = logical
   1

% just call diag() a bazillion times for each channel
timeit(@() advectorize1(A0))

ans = 0.0339

timeit(@() test1(A0))

ans = 0.1056

function B = test1(A0)
H = hankel(1:size(A0,1),size(A0,1):(size(A0,1)+size(A0,2)-1)).';
[~,ind] = sort(H(:));
Y = squeeze(reshape(permute(A0,[2,1,3]),size(A0,1)*size(A0,2),1,[]));
B = Y(ind,:);
end

function B = advectorize1(A)
    A = flipud(permute(A,[2 1 3]));
    s = size(A);
    B = zeros(prod(s(1:2)),s(3));
    for c = 1:s(3)
        Ac = A(:,:,c);
        ipr = 0;
        for k = -(s(1)-1):(s(2)-1)
            d = diag(Ac,k);
            inx = ipr + numel(d);
            B(ipr+1:inx,c) = d;
            ipr = inx;
        end
    end
end

[DGM]

I think that's pretty much as close as anything to what I was shooting for. I'd frustrated myself trying to come up with something using toeplitz() or hankel(), but I'm not exactly clever. It'd be nice if elegant solutions were fast too, but I think at this point I was mostly trying to satisfy that curiosity.

I might play with this for a bit later tonight.

Answer 2

[John D'Errico]

This should not be that difficult. One trick might be to use ndgrid, and then a sort. For example, consider the function adunroll below.

A = reshape(1:12,[3,4]) % input

A = 3×4

1 4 7 10 2 5 8 11 3 6 9 12

adunroll(A)

ans = 12×1

1 4 2 7 5 3 10 8 6 11

function V = adunroll(A)
  % unrolls the general nxm matrix A into a column vector in an anti-diagonal raster form

  % rows and column of A
  [nr,nc] = size(A);
 
  % use ndgrid to turn this into a simple sort problem
  [R,C] = ndgrid(1:nr,1:nc);

  % Look carefully at how M was constructed
  M = [reshape(R+C,[],1),R(:)];

  % sorting the rows of M will result in the order we want
  [~,ind] = sortrows(M);

  % just use those indices to extract the elements of A
  V = A(ind);
end

I won't test the time required, but the above code should be eminently reasonable in terms of time required, mainly becaue it needs only a few very simple function calls. The nice thing is, the code is what I would call elegant, clean and simple to follow. It also has no need to worry about the general shape of A. Square or rectangular in any form will still work. All of that makes the code good, since it will be easy to debug and maintain.

Look carefully at how I created the array M, as that is the key to making this work for any size of array. For example, look at the matrix (R+C). Does that give you a hint at what it was done?

The trick of building on ndgrid is a nice one to remember.

I was going to compare the time required, but I see that your code wants to take an nxmx3 array, but your question was specific to a nxm array. The trick here would be simple enough. Only compute the indices ONCE, then apply them to each panel of the 3 dimensional array.

[DGM]

Thank you for another good solution.

Extending it to operate pagewise is easy enough. It's not the fastest, but it's relatively transparent to read. I think Paul's suggestion is easier to read in some ways, but maybe that's just because I've spent more time playing with similar attempts.

% build 1000x1000x3 test array
N = 1000;
s = [N N];
A0 = reshape(1:prod(s),s);
A0 = repmat(A0,[1 1 3]);

timeit(@() adunroll(A0))

ans = 0.1715

function V = adunroll(A)
    % unrolls the general nxm matrix A into a column vector in an anti-diagonal raster form
    % rows and column of A
    [nr,nc,np] = size(A);
    
    % use ndgrid to turn this into a simple sort problem
    [R,C] = ndgrid(1:nr,1:nc);
    % Look carefully at how M was constructed
    M = [reshape(R+C,[],1),R(:)];
    % sorting the rows of M will result in the order we want
    [~,ind] = sortrows(M);
    
    % just use those indices to extract the elements of A
    %V = A(ind);
    % one output column per vectorized page
    V = A((ind + prod([nr nc])*(0:np-1)));
    V = reshape(V,[],np);
end

Answer 3

[Chris]

I believe this works for tall arrays and square arrays. I haven't yet figured out how to compensate if N>M.

It's not as fast as some of yours, but if you're expecting some standard array sizes you could precalculate idxs, and the final step doesn't take all that long. (Also, I hear timeit() is better for timing things. I was trying tic and toc, but every time I re-ran the code in this editor, the time increased)

N = 1000;
s = [N N];
A0 = reshape(1:prod(s),s);
A0 = repmat(A0,[1 1 3]);

fun1 = @() doCalc(A0);
timeit(fun1)

ans = 0.1012

idxs = returnidxs(A0);
fun2 = @() reshape(A0(idxs(:)),[],3);
timeit(fun2)

ans = 0.0254

function doCalc(A0)
    [R,C,np] = size(A0);
    bignum = max(R,C);
    lilnum = min(R,C);
    mult = R*C;
    
    xmat = ones(bignum).*(1:bignum)';
    r1 = triu(ones(C).*(1:C)');
    r2 = tril(xmat,-1)-tril(xmat,-lilnum-1);
    rowcoords = [r1(:);r2(:)];
    rowcoords = rowcoords(rowcoords~=0);
    
    ymat = (xmat(1:lilnum,:))';
    c1 = flipud((tril(ymat(1:end-1,:)))');
    c2 = flipud(triu(ymat)');
    colcoords = [c1(:);c2(:)];
    colcoords = colcoords(colcoords~=0);
    
    idxs = sub2ind(size(A0),rowcoords,colcoords)+mult*(0:np-1);
    
    B = reshape(A0(idxs(:)),[],np);
end

function idxs = returnidxs(A0)
    [R,C,np] = size(A0);
    bignum = max(R,C);
    lilnum = min(R,C);
    mult = R*C;
    
    xmat = ones(bignum).*(1:bignum)';
    r1 = triu(ones(C).*(1:C)');
    r2 = tril(xmat,-1)-tril(xmat,-lilnum-1);
    rowcoords = [r1(:);r2(:)];
    rowcoords = rowcoords(rowcoords~=0);
    
    ymat = (xmat(1:lilnum,:))';
    c1 = flipud((tril(ymat(1:end-1,:)))');
    c2 = flipud(triu(ymat)');
    colcoords = [c1(:);c2(:)];
    colcoords = colcoords(colcoords~=0);
    
    idxs = sub2ind(size(A0),rowcoords,colcoords)+mult*(0:np-1);
end

[DGM]

I edited the post to use timeit(). That way the tool and machine are at least the same. Your method is pretty much right in the middle of the pack like you said. The output indexing and reshaping eats a lot of time just like in the other cases.

I haven't tried messing around with non-square stuff yet, but returnidxs() tries to return a 3-col matrix regardless of how many channels A has. I haven't wrapped my head around it enough to figure out where that's happening, but that's another thing to consider. I did say (M*N)x3, but I really should've said (M*N)x(numchannels).

[Chris]

Oh, I assumed these were going to be images. I changed two lines to fix the numpages thing.

The indexing of the first page produces a pair of patterns--Row indexes count up until one edge is encountered:

1, (1 2), (1 2 3), .. (1 .. N) -> (2 .. N), (3 .. N), .., N

and column indexes count down

1, (2 1), (3 2 1), .. (N .. 1) -> (N .. 2), (N .. 3), .., N

When it's a rectangle, there's a period in the middle where one index repeats the same series of numbers, corresponding to 1 .. lilnum or lilnum .. 1, while the other index shifts down the line.

This was a late night attempt at capturing that behavior, but in the end it's essentially a bunch of diags (tril, triu) as well.

[DGM]

You were correct to assume that I'm aiming to work with images, but my assumed definition of an image may have 1, 2, 3, 4, or 6 channels.

[Chris]

This is about as well as I can do, and I don't think it's any more elegant than your code or Paul's. Those diags/hankels are fast.

N = 1000;
s = [N N];
A0 = reshape(1:prod(s),s);
A0 = repmat(A0,[1 1 3]);

fun = @() doCalc(A0);
timeit(fun)

ans = 0.0567

function B = doCalc(A0)

    [R,C,np] = size(A0);
    bignum = max(R,C);
    lilnum = min(R,C);
    delta = bignum-lilnum;

    basevec1 = (1:lilnum)';

    starty = triu(repmat(basevec1,1,lilnum));    
    startx = flipud(starty);

    if C>=R

        midx = repmat((lilnum+1:-1:2)',1,delta)+(0:delta-1);
        midy = repmat((1:lilnum)',delta,1);

        basevec2 = (bignum:-1:1+delta)';
        endy = tril(repmat(basevec1,1,lilnum));
        endx = fliplr(triu(repmat(basevec2,1,lilnum)));

    else
        
        midx = repmat((lilnum:-1:1)',1,delta);
        midy = repmat((2:lilnum+1)',1,delta)+(0:delta-1);

        basevec2 = (1+delta:bignum)';
        endy = tril(repmat(basevec2,1,lilnum));
        endx = fliplr(triu(repmat(flipud(basevec1),1,lilnum)));
            
    end

    colidx = [startx(startx~=0);midx(:);endx(endx(:,2:end)~=0)];
    rowidx = [starty(starty~=0);midy(:);endy(endy(:,2:end)~=0)];

    idxs = sub2ind([R,C],rowidx,colidx)+R*C*(0:np-1);

    B = A0(idxs);

end

[DGM]

That's better still. Your method isn't winning awards on the "elegance" front, but it's faster than Paul's method and half of what I came up with.

Yeah. Even before I posted the question, I was kind of baffled that calling diag() repeatedly could be a practical option.

[Chris]

I think diagonals are part of LAPACK, so it would be difficult to find something more efficient.

Shave off another hundredth if you don't need the indexes, by combining the last few lines:

B = A0(sub2ind([R,C],...
    [starty(starty~=0);midy(:);endy(endy(:,2:end)~=0)],...
    [startx(startx~=0);midx(:);endx(endx(:,2:end)~=0)])+R*C*(0:np-1));