# A better way to do antidiagonal matrix vectorization?

6 ビュー (過去 30 日間)
DGM 2021 年 11 月 12 日
コメント済み: DGM 2021 年 11 月 14 日
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
ans = 0.0423
% use an index array and only call diag() one bazillion times
ans = 0.0519
% build index array, nan pad, vector circshift
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
ans = 0.1761
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
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
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
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?

サインインしてコメントする。

### 採用された回答

Paul 2021 年 11 月 12 日

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]);
B2 = test1(A0);
isequal(B1,B2) % check
ans = logical
1
% just call diag() a bazillion times for each channel
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
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 件のコメント表示非表示 なし
DGM 2021 年 11 月 12 日
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.

サインインしてコメントする。

### その他の回答 (2 件)

John D'Errico 2021 年 11 月 13 日

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
ans = 12×1
1 4 2 7 5 3 10 8 6 11
% 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 件のコメント表示非表示 なし
DGM 2021 年 11 月 14 日
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]);
ans = 0.1715
% 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

サインインしてコメントする。

Chris 2021 年 11 月 12 日

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 件のコメント表示非表示 5 件の古いコメント
Chris 2021 年 11 月 13 日
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));

サインインしてコメントする。

R2019b

### Community Treasure Hunt

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

Start Hunting!