Fastest way to sum the upper triangle

Question

I feel like this is an recurring question: if there's a symmetric matrix whose diagonal is not all 0, how could I get the sum of the part of it that's above the diagonal as fast as possible?

Small working example:

a = (# + Transpose[#]) &[RandomReal[{-1, 1}, {1000, 1000}]];

(Total[a, 2] - Tr[a])/2 // AbsoluteTiming (*{0.004000,721.277}*)

Sum[a[[i, j]], {i, 1, 1000}, {j, i + 1, 1000}] // AbsoluteTiming (*{1.254072,721.277}*)

The latter is so much slower than the former. But calculating the sum twice still seem to be a bit unnecessary, could it be faster?

Answer 1

I found this question quite interesting, so I thought I would collect the answers contributed in comments for future reference and to have the question appear as answered in search. I generated a slightly bigger matrix to play with, and minimally modified the code to render it independent of the size of the matrix. I also compared timings of each method to have a homogeneous set at least on my system (MMA 10.2 on Win7-64bit).

SeedRandom[23]
a = (# + Transpose[#])& [RandomReal[{-1, 1}, {5000, 5000}]];

---

If the matrix is symmetric, as was the case in the OP and in the example above, then taking advantage of its structure yields the fastest solution so far, already originally proposed by arax in the question:

(Total[a, 2] - Tr[a])/2 // RepeatedTiming
(* Out: {0.025, -1769.6} *)

---

In the general case of a non-symmetrical matrix, however, the above solution of course wouldn't work. The Toad proposed a very compact solution using UpperTriangularize. Timings show that Flatten is bad for speed; much more expedient is to use a levelspec option for Total.

Total@Flatten@UpperTriangularize[a, 1] // RepeatedTiming
Flatten@a; // RepeatedTiming

Total[UpperTriangularize[a, 1], 2] // RepeatedTiming

(* Out: 
{0.338, -1769.6}
{0.20, Null}

{0.134, -1769.6}
*)

Similarly one could use LowerTriangularize to get the sum of the lower triangular portion of a non-symmetric matrix: Total[LowerTriangularize[a, -1], 2].

---

Kale then proposed a faster approach based on explicit part extraction. The version proposed in the comment actually sums the lower triangular members of the original matrix, which made no difference in the case of a symmetric matrix, but of course would matter in the general case. Two versions are presented here, for upper and lower triangles (which also required a slight adjustment to the Range used):

(* Upper triangular version *)
Total[a[[#, # + 1 ;; -1]] & /@ Range[1, Length[a]], 2] // RepeatedTiming
(* Out: {0.0940, -1769.6} *)

(* Lower triangular version *)
Total[a[[#, 1 ;; # - 1]] & /@ Range[1, Length[a]], 2] // RepeatedTiming
(* Out: {0.104, -1769.6} *)

Although faster in this case, the intermediate list of parts generated in this approach is ragged, which may lead to unpacking of packed arrays (see e.g. Leonid's discussion on the topic).

---

N.J. Evans proposed a simple modification to Kale's solution that yields a significantly faster method, by simply applying Total twice in a row, rather than using a levelspec:

Total[Total@a[[#, # + 1 ;; -1]] & /@ Range[1, Length[a]]] // RepeatedTiming
(* Out: {0.04, -1769.6} *)

This works exceedingly well for the part-extraction approach, where it may have to do with the fact that the outer Total is presented with the much simpler job of summing over a list, thanks to the inner Total, which may also have had the advantage of quick mapping.

In support of this hypothesis is also the fact that the same approach does not yield faster results in the *Triangularize approach.

---

Finally, the approach based on Sum presented in the OP as a slow comparison is, indeed, very slow:

Sum[a[[i, j]], {i, 1, Length[a]}, {j, i + 1, Length[a]}] // RepeatedTiming
(* Out: {18.3, -1769.6} *)

Answer 2

For a non-symmetric real matrix you can consider using LibraryLink to speed things up. It still won't be as fast as the Total/Tr answer, but it may be useful otherwise (call this C program SumUpperTriangle.c):

#include "WolframLibrary.h"

DLLEXPORT int SumUpperTriangle(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res) {

 /* Variable declarations */

 int err = LIBRARY_NO_ERROR;

 MTensor in;
 mint const *in_dims;
 mreal *in_data;

 mreal sum;
 mint rows, row, cols, col;

 /* Code section */

 /* Get the input matrix */
 in = MArgument_getMTensor(Args[0]);
 in_dims = libData->MTensor_getDimensions(in);
 in_data = libData->MTensor_getRealData(in);

 rows = in_dims[0];
 cols = in_dims[1];
 sum = 0.0;
 for(row=0; row<rows; row++) {
  for(col=row+1; col<cols; col++) {
   sum += in_data[col+cols*row];
  }
 }

 MArgument_setReal(Res,sum);
 return err;

}

To compile this C code and use SumUpperTriangle as a Wolfram language function you can run the following code (call this code SumUpperTriangle.wl):

Needs["CCompilerDriver`"];
lib = CreateLibrary[{"SumUpperTriangle.c"}, "SumUpperTriangle" ];
SumUpperTriangle = LibraryFunctionLoad[lib, "SumUpperTriangle", {{Real, 2}}, Real ];

Now we can load the code above:

Get["SumUpperTriangle.wl"]

And compare it to the Total/Tr solution:

aa := (# + Transpose[#]) &[RandomReal[{-1, 1}, {1000, 1000}]];
a = aa;
SumUpperTriangle[a] == (Total[a, 2] - Tr[a])/2
(* Gives: True *)

Here I am using := to define aa so that aa is different every time you call it. This is useful to avoid timings that include numerical caching results.

Now let's compare timings with the LibraryLink function:

Total[Table[a = aa; First[AbsoluteTiming[SumUpperTriangle[a]]], {100}]] 
(* gives 0.32099 seconds on my machine *)

And the original solution:

Total[Table[a = aa; First[AbsoluteTiming[(Total[a, 2] - Tr[a])/2]], {100}]]
(* gives 0.102354 seconds on my machine *)

So the Total/Tr solution is still faster, but the SumUpperTriangle will work for any real matrix, symmetric or not. And SumUpperTriangle is of course faster than the Sum[...] approach...

For more information on LibraryLink, see the user guide or these introductory examples on the Wolfram community.