A better way to do antidiagonal matrix vectorization?

12 views (last 30 days)
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?

Accepted Answer

Paul
Paul on 12 Nov 2021
Edited: Paul on 12 Nov 2021
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
  1 Comment
DGM
DGM on 12 Nov 2021
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.

Sign in to comment.

More Answers (2)

John D'Errico
John D'Errico on 13 Nov 2021
Edited: John D'Errico on 13 Nov 2021
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.
  1 Comment
DGM
DGM on 14 Nov 2021
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

Sign in to comment.


Chris
Chris on 12 Nov 2021
Edited: Chris on 12 Nov 2021
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
  6 Comments
DGM
DGM on 13 Nov 2021
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
Chris on 13 Nov 2021
Edited: Chris on 13 Nov 2021
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));

Sign in to comment.

Categories

Find more on Loops and Conditional Statements in Help Center and File Exchange

Products


Release

R2019b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!