# Take upper triangular part of matrix

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 = DeveloperFromPackedArray[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 '18 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 '18 at 19:35

## 2 Answers

Is using an internal, undocumented symbol acceptable?

r1 = StatisticsLibraryUpperTriangularMatrixToVector[mat]; //RepeatedTiming
r2 = StatisticsLibraryUpperTriangularMatrixToVector[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 '18 at 15:12
• 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 '18 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 '18 at 16:51
• It seems that for Mma 10.0, this symbol no longer exists. – Denis Cousineau Mar 4 at 20:57

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 '18 at 14:11