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.

Answer 1

[Andrei Bobrov]

"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(:));

[Jon Paul Janet]

This is the most elegant answer I think. Unfortunately it is pretty similar in run time to what I am currently doing, but it's certainly neater.

[Walter Roberson]

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]

thank you very much sir.... very elegant solution

[Hongbo Zhao]

Very elegant indeed.

Answer 2

[Walter Roberson]

[r, c] = size(A);
sum(A(r : r-1 : end))

[Stephen23]

+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]

For non-square matrices,

sum( X(r:r-1:r*(r-1)+1) )

Answer 3

[Thorsten]

You could try

   S = sum(sum(A - diag(diag(A))));

[Jon Paul Janet]

I'm sorry if I was not clear, but the anti-diagonal elements are those on the NE-SW diagonal, like this, except generalized in the same way as diag(A,1) to diag(A,N/2). In other words, all of the elements on any NE-SW diagonal

Answer 4

[Sean de Wolski]

diag(flip(A))

[Walter Roberson]

I don't think this is going to be the most efficient...