Pfthb

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.

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

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