File Exchange

## sumDiag

version 1.4 (1.59 KB) by

Sum over each diagonal (or anti-diagonal) in a matrix

Updated

Efficient and compact code for summing each diagonal (or anti-diagonal) in a matrix without using a for-loop.
Works well for large matrices.
For 3D matrix input A, the sum of diagonals of A(:,:,k) are returned in
sumMat(:,k). The script is typically faster than a for-loop based approach when A is 3D.
The code is most efficient for wide or tall matrices. Inline the code when it is used as a part of an iterative algorithm, to avoid recomputing constant indexing matrices.

Note that a for-loop implementation with the diag() function can be faster, and has a lower memory requirement, especially in the square 2D case.

Ton D

### Ton D (view profile)

For 2D matrices if you want the anti-diagonal sums the following code is much faster (and not hard to adapt for for ND matrices or diagonal sums):

m=size(A,1);
n=size(A,2);

tmp = [A;zeros(n)];
tmp = tmp(:);
res = sum(reshape(tmp(1:end-n),m+n-1,n),2);

Marcus Björk

### Marcus Björk (view profile)

That's true Sven. I added 3D matrix support now, if anyone has an application for that.
In my case the code was used in an iterative algorithm, in which case the matrix is not known beforehand and cannot be stored in a 3D matrix. I figured the code was nice enough to make a function out of it and upload it here.

Sven

### Sven (view profile)

Ha, but if you used the sumDiags() function itself multiple times you wouldn't actually get your speedup :)

I'd say if that's the goal then perhaps you can rewrite sumDiags to take in a 3D matrix and give the output as the sum of diags of each sheet.

In the little tests I ran I found that the sum(diag()) loop was about twice as fast for square matrices and about equal for really wide (10x10000) matrices.

Oh, and yep you're right that rot90() is the way to go for anti-diagonals :)

Marcus Björk

### Marcus Björk (view profile)

Forgot to mention: If you need to compute the sum of diagonals for several matrices of the same size (in another loop), which was the application at hand, this for-loop free implementation is also faster (well, the inlined version). See below:

Test code:
%Input
K=1000;
A=randn(1000,100);

%Common
[N,M]=size(A);

%% Method 1
tic;
d = 1-N:M-1;
out = zeros(M+N-1,1,class(A));

for k=1:K
for p = 1:length(d)
out(p) = sum(diag(A,d(p)));
end
A=randn(N,M);
end
toc

%% Method 2
tic;
Amod=zeros(N+M-1,M);
logVec = [false(M-1,1);true(N,1);false(M-1,1)];
indMat = bsxfun(@plus, (1:M+N-1)',0:M-1);
logMat = logVec(indMat);

for k=1:K
Amod(logMat)=A;
sumVec=sum(Amod,2);
A=randn(N,M);
end
toc

>>
Elapsed time is 2.632392 seconds.
Elapsed time is 1.307457 seconds.

Marcus Björk

### Marcus Björk (view profile)

I had the for-loop implementation already, which is trivial, but was posed with the problem of solving it without a for-loop. Hence, I wrote this code. (Which could perhaps be optimized further?)
Didn't actually check which was faster, and I should probably make a note of this in the description.

Furthermore, you don't get the anti-diagonals by transposing. But rot90(A) works and is almost as fast as using the extra input.

Sven

### Sven (view profile)

Hi Marcus,

Nice entry, but sometimes a for-loop can be useful. The function below is simpler, uses less memory (no extra matrices or masks needed), and is actually more efficient (faster).

function out = sumDiag(A)
d = 1-size(A,1):size(A,2)-1;
out = zeros(size(d),class(A));
for k = 1:length(d)
out(k) = sum(diag(A,d(k)));
end

I chose the opposite convention, whereby the first element of the output is the lower left rather than upper right. I think this fits in with MATLAB's convention of numbering diagonals as lower to the left and higher to the right, but you can change it easily if you'd like on the first line.

Also, I think that the anti-diagonals case should really just be up to the user who can call sumDiag(A) or sumDiag(A'), as that is simpler for the same result than (possibly) providing an extra argument.