Code covered by the BSD License  

Highlights from
sumDiag

Be the first to rate this file! 13 Downloads (last 30 days) File Size: 1.59 KB File ID: #44162

sumDiag

by

 

01 Nov 2013 (Updated )

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

| Watch this File

File Information
Description

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.

Required Products MATLAB
MATLAB release MATLAB 8.2 (R2013b)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (5)
05 Nov 2013 Marcus Björk

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.

03 Nov 2013 Sven

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 :)

02 Nov 2013 Marcus Björk

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.

02 Nov 2013 Marcus Björk

Thanks for your comment Sven!
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.

02 Nov 2013 Sven

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.

Updates
04 Nov 2013

Added comment on the speed of the implementation compared to using a for-loop approach.

05 Nov 2013

Added support for 3D matrices as input.

Contact us