11
$\begingroup$

I would like to extract the upper triangular part of a square matrix into a flat list. I am looking for fast ways to do this. I am primarily interested in solutions that are fast for non-packed arrays as well. For packed arrays, a compiled solution will always be very fast.

Here are three successively faster implementations to get the ball rolling.

takeUpper1[mat_?SquareMatrixQ] := 
  Join @@ Table[mat[[i, j]], {i, Length[mat]}, {j, i + 1, Length[mat]}];

takeUpper2[mat_?SquareMatrixQ] := 
  Extract[mat, Subsets[Range@Length[mat], {2}]]

takeUpper3[mat_?SquareMatrixQ] := 
  Join @@ Pick[mat, UpperTriangularize[ConstantArray[1, Dimensions[mat]], 1], 1]

(* added later, also slower than takeUpper3 *)
takeUpper4[mat_?SquareMatrixQ] := 
  Join @@ Table[mat[[i, i + 1 ;;]], {i, Length[mat]}]

Benchmarks:

mat = RandomReal[1, {1000, 1000}];    
mat2 = Developer`FromPackedArray[mat];

Packed:

res1 = takeUpper1[mat]; // RepeatedTiming
(* {0.218, Null} *)

res2 = takeUpper2[mat]; // RepeatedTiming
(* {0.0840, Null} *)

res3 = takeUpper3[mat]; // RepeatedTiming
(* {0.017, Null} *)

res1 == res2 == res3
(* True *)

Non-packed:

res1 = takeUpper1[mat2]; // RepeatedTiming
(* {0.839, Null} *)

res2 = takeUpper2[mat2]; // RepeatedTiming
(* {0.0949, Null} *)

res3 = takeUpper3[mat2]; // RepeatedTiming
(* {0.018, Null} *)

res1 == res2 == res3
(* True *)

You are also welcome to suggest intuitive names for such a function. The final function will work also on non-square matrices and will have a second argument similar to that of UpperTriangularize.

$\endgroup$
2
  • $\begingroup$ Are you interested in only positive second arguments? $\endgroup$
    – Carl Woll
    Nov 9, 2018 at 17:08
  • $\begingroup$ @CarlWoll Also negative ones, but I think that would be too much for this question. This question was only for the square case, no second argument, and focusing on performance. I was useful to learn about UpperTriangularMatrixToVector. $\endgroup$
    – Szabolcs
    Nov 11, 2018 at 19:35

2 Answers 2

14
$\begingroup$

Is using an internal, undocumented symbol acceptable?

r1 = Statistics`Library`UpperTriangularMatrixToVector[mat]; //RepeatedTiming
r2 = Statistics`Library`UpperTriangularMatrixToVector[mat2]; //RepeatedTiming
r3 = takeUpper3[mat]; //RepeatedTiming
r4 = takeUpper3[mat2]; //RepeatedTiming

r1 === r2 === r3 === r4

{0.00039, Null}

{0.0027, Null}

{0.017, Null}

{0.018, Null}

True

$\endgroup$
4
  • $\begingroup$ Wow! That's fast. Have to remember that one. $\endgroup$ Nov 9, 2018 at 15:12
  • 3
    $\begingroup$ Yes, it is acceptable. But it obsoletes all the work I did on this so far ... This stuff should really be documented. $\endgroup$
    – Szabolcs
    Nov 9, 2018 at 15:30
  • $\begingroup$ Just a note that this function does not simply re-pack mat2 and use the implementation for packed arrays. I tested this by trying it on matrices of strings. It is still fast. $\endgroup$
    – Szabolcs
    Nov 9, 2018 at 16:51
  • $\begingroup$ It seems that for Mma 10.0, this symbol no longer exists. $\endgroup$ Mar 4, 2021 at 20:57
5
$\begingroup$

I tried the following; for packed arrays, it appears to be on par with takeUpper3, but it needs twice the time for unpacked arrays. So I think, your trick using Pick is already pretty good.

LinearToTriangularIndexing[k_?VectorQ, n_Integer] := Module[{i, j},
   i = n - 1 - Floor[Sqrt[4. n (n - 1) - 8. k + 1.]/2.0 - 0.5];
   j = Subtract[
     k + i + Quotient[Subtract[n + 1, i] Subtract[n, i], 2], 
     Quotient[n (n - 1), 2]];
   Transpose[{i, j}]
   ];

takeUpper5[mat_?SquareMatrixQ] := With[{n = Length[mat]},
  Extract[mat, LinearToTriangularIndexing[Range[1, n (n - 1)/2], n]]
  ]
$\endgroup$
1
  • $\begingroup$ For packed arrays, I have a C solution that runs un 0.0007 s for this matrix. This is why I am focused on things that are not packable. $\endgroup$
    – Szabolcs
    Nov 9, 2018 at 14:11

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.