Most efficient way to sum anti-diagonal elements
显示 更早的评论
Inside an ODE that I am solving (ode15s), I need to do the following. Let the state vector by y (Nx1), and for some fixed, sparse symmetric A (NxN), I need to have that
y' = [sum of anti-diagonal elements of (diag(y)*A*diag(y))] + f(t)
for some forcing f. The problem is, for large N (10k + ) this is pretty slow as the solver takes many steps. Currently I have calculated a logical index mask (outside the ode) that gives me the elements I want, and my code (inside the d.e.) is:
A = spdiags(y,0,N,N)*A*spdiags(y,0,N,N); %overwrite A in-place with y*A*y
working_matrix=sparse([],[],[],2*N-1,N,N*ceil(N/2)); % sparse allocation
working_matrix(index_map)=A; %strip antidiagonals to columns
this gives me working_matrix, which has the antidiagonal elements of (diag(y)*A*diag(y) which I can just sum over. However, 99% of my runtime is spent on the line "working_matrix(index_map)=A". Any speedup on this line would save me a lot of time. "index_map" is a (2*N-1)xN logical array that pulls out the correct elements, based on this work here.
Is there a better way? I can pull of the antidiagonal elements of A outside the solver and pass the matrix that has the antidiagonal elements as rows, but then I still need the same construction applied to y*y' to get the matching elements of y - unless there is a better way to do this?
I am running R2015b if that matters.
采纳的回答
Andrei Bobrov
2016-7-12
"all of the elements on any NE-SW diagonal"
a = randi(50,6);
[m,n] = size(a);
idx = hankel(1:m,m:(n-1)+m);
out = accumarray(idx(:),a(:));
4 个评论
Jon Paul Janet
2016-7-12
Walter Roberson
2016-7-12
This is code for working on all of the diagonals at the same time. If you do not need all of them, then this will be much less efficient than my solution.
Salvatore Russo
2018-5-3
thank you very much sir.... very elegant solution
Hongbo Zhao
2018-10-10
Very elegant indeed.
更多回答（3 个）
Walter Roberson
2016-7-11
[r, c] = size(A);
sum(A(r : r-1 : end))
2 个评论
+1 for a very neat solution. Note the indices must end with end-1:
>> X = [1,2,3;4,5,6;7,8,9]
X =
1 2 3
4 5 6
7 8 9
>> [r,c] = size(X);
>> X(r:r-1:end-1)
ans =
7 5 3
>> X(r:r-1:end)
ans =
7 5 3 9
The code should be:
>> sum(X(r:r-1:end-1))
ans = 15
Also note that this will only work for square matrices (and a few other cases) but not for any general rectangular matrix.
Walter Roberson
2016-7-12
For non-square matrices,
sum( X(r:r-1:r*(r-1)+1) )
Sean de Wolski
2016-7-11
diag(flip(A))
1 个评论
Walter Roberson
2016-7-12
I don't think this is going to be the most efficient...
类别
在 帮助中心 和 File Exchange 中查找有关 Operating on Diagonal Matrices 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
