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.

  • Are you interested in only positive second arguments? – Carl Woll Nov 9 at 17:08
  • @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. – Szabolcs Nov 11 at 19:35
up vote 11 down vote accepted

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

  • Wow! That's fast. Have to remember that one. – Henrik Schumacher Nov 9 at 15:12
  • 2
    Yes, it is acceptable. But it obsoletes all the work I did on this so far ... This stuff should really be documented. – Szabolcs Nov 9 at 15:30
  • 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. – Szabolcs Nov 9 at 16:51

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]]
  ]
  • 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. – Szabolcs Nov 9 at 14:11

Your Answer

 
discard

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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