Take upper triangular part of matrix

Question

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.

Answer 1

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

Answer 2

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