14
$\begingroup$

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?

$\endgroup$
  • 8
    $\begingroup$ Another way Total@Flatten@UpperTriangularize[a, 1] (or Total[..., 2] w/o Flatten). Probably not as fast as the Total - Tr route, but easier to read perhaps. $\endgroup$ – rm -rf Aug 6 '15 at 23:03
  • $\begingroup$ @TheToad: Beat me to it, you clever frog...+1 $\endgroup$ – ciao Aug 6 '15 at 23:04
  • $\begingroup$ @TheToad Great! (Un?)surprisingly Flatten took half of the time in the first expression. Both are fast enough to sum a matrix the size of my available RAM. $\endgroup$ – arax Aug 6 '15 at 23:07
  • $\begingroup$ Total[a[[#, ;; # - 1]] & /@ Range[2, 1000], 2] // AbsoluteTiming isn't too bad. $\endgroup$ – kale Aug 7 '15 at 0:14
  • $\begingroup$ @arax If one of the two answers you received seems to appropriately answer your question, you might want to consider accepting it for future reference. $\endgroup$ – MarcoB Aug 18 '15 at 21:51
8
$\begingroup$

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} *)
$\endgroup$
  • 1
    $\begingroup$ A slight modification to Kale's answer gives me about double the speed: Total[Total@a[[#, # + 1 ;; -1]] & /@ Range[1, Length[a]]] // AbsoluteTiming $\endgroup$ – N.J.Evans Aug 10 '15 at 20:58
  • $\begingroup$ @N.J.Evans That's remarkable! I had tried a Total@Total version on the *Triangularize method, and in that case it's slower than Total[#,2]. I hope you won't mind if I add your comment to the answer. $\endgroup$ – MarcoB Aug 10 '15 at 21:07
  • $\begingroup$ Not at all. I was surprised too, but it makes sense. The innermost Total reduces the sub-arrays to single numbers during the mapping, and then the outer Total sums over a 1-D list. $\endgroup$ – N.J.Evans Aug 10 '15 at 21:11
  • 1
    $\begingroup$ @N.J.Evans Added. Thanks for bringing this up! It's quite the improvement speed-wise. I also commented on the possible reason for the speedup. Feel free to edit if you don't think it appropriate. $\endgroup$ – MarcoB Aug 10 '15 at 21:19
7
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Very nice; I had wondered how much this approach would yield. It is worth noting though that Total uses compensated summation (the Kahan algorithm), which isn't used here. Therefore, the roundoff error can be larger. However, Kahan's algorithm is not difficult to implement, so your point stands. In any case, it may be useful to unroll the inner loop a little, and perhaps to accumulate the row sums in parallel using OpenMP. $\endgroup$ – Oleksandr R. Aug 10 '15 at 22:18
  • $\begingroup$ Hi Arnoud, you've posted many LibraryLink examples. Could you take a look at this question? Unrelated: you might also be interested in this. $\endgroup$ – Szabolcs Aug 19 '15 at 20:21
  • $\begingroup$ LTemplate would reduce the amount of code one needs to write like this. Unfortunately on my machine the LibraryLink function is twice slower than the builtin approach. But that's not due to LTemplate, it's the same if I use your code. $\endgroup$ – Szabolcs Aug 19 '15 at 20:37
  • 2
    $\begingroup$ I never realized that WolframLibrary_initialize, etc. were not mandatory in every library. $\endgroup$ – Szabolcs Aug 19 '15 at 20:52
  • $\begingroup$ I think it is good practice to have them, especially in production code. But for simple examples like these, everything seems to work fine without them. $\endgroup$ – Arnoud Buzing Aug 20 '15 at 15:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.