# Populate an upper triangular matrix from a vector of elements

What is the best way to populate an upper triangular (or alternatively: lower triangular or symmetric) matrix from a vector of elements?

This example input:

{1, 2, 3, 4, 5, 6}


should give

{{0, 1, 2, 3},
{0, 0, 4, 5},
{0, 0, 0, 6},
{0, 0, 0, 0}}


or alternatively

{{0, 1, 2, 3},
{1, 0, 4, 5},
{2, 4, 0, 6},
{3, 5, 6, 0}}


or even

{{0, 0, 0, 0},
{1, 0, 0, 0},
{2, 3, 0, 0},
{4, 5, 6, 0}}


This problem has many solutions of course. I'm looking for the best (criteria: most elegant/readable, shortest, fastest) ones.

The input has length $\binom{n}{2} = \frac{n(n-1)}{2}$. Consider $n$ known. Feel free to use any version 10 functionality.

• @MrWizard I hope you don't mind, I know we went through this once in the past. Perhaps v10 has new ways. Note: I do have a number of solutions to the problem (including MrW's), but I chose not to share them in order not to bias the answers. Commented Jul 23, 2014 at 14:04
• The elements are machine integers only? Or...? Commented Jul 23, 2014 at 14:14
• @OleksandrR. I was thinking of a general solution, but actually I will only use machine integers in practice. A fast solution using machine integers only will definitely be useful. Commented Jul 23, 2014 at 14:16

Here is a semi-imperative function to create an upper-triangular array:

upperTriangular[v_] := upperTriangular[v, (1 + Sqrt[1 + 8*Length@v])/2]

upperTriangular[v_, n_] := Module[{i = 0}, Array[If[# >= #2, 0, v[[++i]]]&, {n, n}]]


The function is expressed using two definitions for a reason that will become clear in a moment. Here it is in action:

upperTriangular @ Range @ 6

(* {{0,1,2,3},{0,0,4,5},{0,0,0,6},{0,0,0,0}} *)


If we know the square matrix size n ahead of time, and we intend to perform many such triangularizations, then it might be useful to precompile the transformation function. The following function does that, using the two-argument form of upperTriangular to write the matrix construction code:

upperTriangulizer[n_Integer, type_:_Integer] :=
Module[{v, r}
, Quiet[With[{r = upperTriangular[v, n]}, Compile[{{v, type, 1}}, r]], Part::partd]
]


Edit see the note below if this expression produces a warning

Here it is in action:

upperTriangulizer[5] @ Range[10]


(* {{0,1,2,3,4},{0,0,5,6,7},{0,0,0,8,9},{0,0,0,0,10},{0,0,0,0,0}} *)

The result is a packed array:

DeveloperPackedArrayQ @ %

(* True *)


We can control the data type of the matrix elements:

upperTriangulizer[4, _Real] @ Range[6]

(* {{0.,1.,2.,3.},{0.,0.,4.,5.},{0.,0.,0.,6.},{0.,0.,0.,0.}} *)


An interesting feature of the generated function is that it involves no loops -- the resultant matrix is created by direct construction. The loops have been completely unrolled, a desirable situation for some applications:

CompiledFunctionToolsCompilePrint@upperTriangulizer[3, _Real]

(*
1 argument
4 Integer registers
9 Real registers
2 Tensor registers
Underflow checking off
Overflow checking off
Integer overflow checking on
RuntimeAttributes -> {}

T(R1)0 = A1
I0 = 0
I2 = 2
I1 = 1
I3 = 3
Result = T(R2)1

1   R0 = Part[ T(R1)0, I1]
2   R1 = Part[ T(R1)0, I2]
3   R2 = Part[ T(R1)0, I3]
4   R3 = I0
5   R4 = I0
6   R5 = I0
7   R6 = I0
8   R7 = I0
9   R8 = I0
10  T(R2)1 = {{R3, R0, R1}, {R4, R5, R2}, {R6, R7, R8}}
11  Return
*)


Note about the warning Optional::opdef

Mathematica versions 9 and earlier issue a warning when the compiled version of upperTriangulizer is defined:

Optional::opdef: "The default value for the optional argument type_:_Integer contains a pattern."

This message is defensive in nature, warning that the construction type_:_Integer is unusual. So unusual, in fact, it is frequently more likely to be a typing error rather than intentional. But in this case, the construction is intentional. We are defining a function that expects a type pattern as an argument, and we are assigning _Integer as the default pattern.

For whatever reason, version 10 no longer issues this warning message. Perhaps WRI is finding such constructions to be more common now? In any event, the message is safe to suppress on earlier versions:

Quiet[
upperTriangulizer[n_Integer, type_:_Integer] :=
Module[{v, r}
, Quiet[With[{r = upperTriangular[v, n]}, Compile[{{v, type, 1}}, r]], Part::partd]
]
, Optional::opdef
]

• Why function upperTriangulizer cannot be executed in my Mathematica V8 enviroment correctly? It gives the warning information:Optional::opdef: "The default value for the optional argument type_:_Integer contains a pattern."
– xyz
Commented Sep 13, 2014 at 9:32
• @Tangshutao I have added a note concerning this message to the bottom of my response. Commented Sep 13, 2014 at 12:07
• I think you can do like upperTriangular[v_, n_] /; Length@v == (n (n - 1))/2:=... to omit the definition upperTriangular[v_] := upperTriangular[v, (1 + Sqrt[1 + 8*Length@v])/2]
– xyz
Commented Sep 14, 2014 at 6:41
• @Tangshutao If I did that, it would place the burden of computing n on the caller and make upperTriangular[v, n] go inert if they got that calculation wrong. If we wanted only one definition, I would keep upperTriangular[n_] and merge the second into it by computing n within the Module. But the separate definitions support both the original use case as well as the precomputation of the transformation function by upperTriangulizer. Commented Sep 14, 2014 at 13:31

Probably not fast, but very simple:

upperTriangular[elements_, n_] :=
SparseArray[
SymmetrizedIndependentComponents[{n, n}, Antisymmetric[{1, 2}]] ->
elements]]

• SymmetrizedIndependentComponents? These function names are getting out of hand. +1 nevertheless as I don't recall seeing that function used before. Commented Jul 23, 2014 at 17:11
• Great! I didn't mention it but I was hoping for something using this functionality. Is StructuredArray useful here? Commented Jul 23, 2014 at 17:25

Update: Recently bumped into StatisticsLibraryVectorToUpperTriangularMatrix which, as the name suggests, seems to do just what is needed. The second argument controls the diagonal elements. The third argument seems to control the matrix dimension; I could get it to work only with appropriate integer based on the length of the first argument.

vecToUTM = StatisticsLibraryVectorToUpperTriangularMatrix[#, #2,
(1 + Sqrt[1 + 8 Length@#])/2]&;

vecToUTM[{1, 2, 3, 4, 5, 6}, 0] // TeXForm


$\left( \begin{array}{cccc} 0 & 1 & 2 & 3 \\ 0 & 0 & 4 & 5 \\ 0 & 0 & 0 & 6 \\ 0 & 0 & 0 & 0 \\ \end{array} \right)$

vecToUTM[{1, 2, 3, 4, 5, 6}, x] // TeXForm


$\left( \begin{array}{cccc} x & 1 & 2 & 3 \\ 0 & x & 4 & 5 \\ 0 & 0 & x & 6 \\ 0 & 0 & 0 & x \\ \end{array} \right)$

vecToUTM[{1, 2, 3, 4, 5, 6} , {a, b, c, d}]


$\left( \begin{array}{cccc} a & 1 & 2 & 3 \\ 0 & b & 4 & 5 \\ 0 & 0 & c & 6 \\ 0 & 0 & 0 & d \\ \end{array} \right)$

ClearAll[f1];
f1 = RotateLeft[PadLeft[InternalPartitionRagged[Range[# (#-1)/2],
Range[#-1, 1, -1]], {#, #}]] &;

Row[MatrixForm /@ {m1 = f1[4], Transpose[m1], m1 + Transpose[m1]}]


ClearAll[f2];
f2 = With[{m = Solve[{a (a - 1)/2 == Length@#, a > 0}, a, Integers][[1, 1, 2]]},
RotateLeft[PadLeft[InternalPartitionRagged[#, Reverse@Range[m - 1]], {m, m}]]] &;

Row[MatrixForm/@ {m1 = f2[Range[6]], Transpose[m1], m1 + Transpose[m1]}]


lst2 = CharacterRange["a", "z"][[;;10]];
Row[MatrixForm /@ {m1 = f2[lst2], Transpose[m1], m1 + Transpose[m1]}]


• Your solution is very efficient
– xyz
Commented Sep 13, 2014 at 14:47
• @Tangshutao, it is - to my surprise- quite fast; @Chip Hurst's compiled Matrixify is twice as fast though. Just crashed Windows trying to do timings for the posted methods:)
– kglr
Commented Sep 13, 2014 at 15:13
• Dear kguler, How can I understand Internal-PartitionRagged, I cannot find the usage of this function in Wolfram Help Documentation
– xyz
Commented Sep 14, 2014 at 6:21
• @Tangshuato, ??InternalPartitionRagged_ gives some info about the function. To see the code you can use _ClearAttributes[InternalPartitionRagged, ReadProtected] followed by ??InternalPartitionRagged.
– kglr
Commented Sep 14, 2014 at 6:27

Here's a procedural, hence compilable solution:

Matrixify = Compile[{{lis, _Integer, 1}},
Module[{n, mat, pl = 1, pm = 2, i = 1},

n = Round[(Sqrt[8 Length[lis] + 1] + 1)/2];
mat = Table[0, {n^2}];

Do[
mat[[pm ;; pm + k - 1]] = lis[[pl ;; pl + k - 1]];
i++;
pl += k;
pm += k + i,
{k, n - 1, 1, -1}
];

Partition[mat, n]

],
CompilationTarget -> "C",
Parallelization -> True,
RuntimeOptions -> "Speed"
];


Here are some results:

Matrixify[Range[6]]
(* {{0, 1, 2, 3}, {0, 0, 4, 5}, {0, 0, 0, 6}, {0, 0, 0, 0}} *)

Matrixify[Range[Binomial[1000, 2]]]; // AbsoluteTiming
(* {0.032628, Null} *)


I realized this is possible to implement rather cleanly using partitionBy as defined here.
This method produces the "or even" form shown as the third example.

in = Range[6];

PadRight @ partitionBy[in, # &] ~ArrayPad~ {{1, 0}, {0, 1}}


$\left( \begin{array}{cccc} 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 2 & 3 & 0 & 0 \\ 4 & 5 & 6 & 0 \\ \end{array} \right)$

## Old answer from Code Review

(Source)

Here are a few methods for your consideration and feedback.

All use:

s = 5
elems = Range[s(s-1)/2]


This one uses the core of my "Dynamic Partition" function. It is the fastest method I know for this problem. Also, perhaps refactoring the code like this makes it more intelligible.

dynP[l_, p_] :=
MapThread[l[[# ;; #2]] &, {{0}~Join~Most@# + 1, #} &@Accumulate@p]

#~PadRight~s & /@ {{}} ~Join~ dynP[elems, Range[s-1]]


This one is slightly shorter, using a different way to construct the indices, but also slightly slower, and perhaps less transparent.

#~PadRight~s & /@
Prepend[
elems~Take~# & /@
Array[{1 + # (# - 1)/2, # (# + 1)/2} &, s - 1],
{}
]


This one is terse, but slow. Legibility is debatable.

SparseArray[Join @@ Table[{i, j}, {i, s}, {j, i - 1}] -> elems, s] // MatrixForm


Here is one I just came up with, and I am quite pleased with it. I think it may be more understandable than most of the ones above.

Take[FoldList[RotateLeft, elems, Range[0, s - 1]] ~LowerTriangularize~ -1, All, s + 1]

s = 4;
lis = Range[s (s - 1)/2];
PadLeft[lis[[1 + # s - (# (1 + #))/2 ;; (1 + #) s - ((# + 1) (2 + #))/
2]], s, 0] & /@ Range[0, s - 1] // MatrixForm

• I think cannot apply to the list of {a,b,c,d,e,f}, it seems that your solution just suitable for the list of Range[n]
– xyz
Commented Sep 14, 2014 at 6:01

## My trial:

Define the function:

upperTriangular[expr_?VectorQ, n_Integer] /; Length@expr == (n (n - 1))/2 :=
Join[Table[0, {n - Length@#}], #] & /@ (Append[
Complement[#2, #1] & @@@
Partition[
Prepend[Take[expr, #] & /@ Accumulate@Reverse@Range[n - 1], {}],2, 1], {}])


Test

expr1=CharacterRange["a", "u"];
upperTriangular[expr, 7] // MatrixForm


$\left( \begin{array}{ccccccc} 0 & \text{a} & \text{b} & \text{c} & \text{d} & \text{e} & \text{f} \\ 0 & 0 & \text{g} & \text{h} & \text{i} & \text{j} & \text{k} \\ 0 & 0 & 0 & \text{l} & \text{m} & \text{n} & \text{o} \\ 0 & 0 & 0 & 0 & \text{p} & \text{q} & \text{r} \\ 0 & 0 & 0 & 0 & 0 & \text{s} & \text{t} \\ 0 & 0 & 0 & 0 & 0 & 0 & \text{u} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \end{array} \right)$

The test of OP

expr2=Range@6;
upperTriangular[expr2, 4] // MatrixForm


$\left( \begin{array}{cccc} 0 & 1 & 2 & 3 \\ 0 & 0 & 4 & 5 \\ 0 & 0 & 0 & 6 \\ 0 & 0 & 0 & 0 \\ \end{array} \right)$

### Testing the efficiency

upperTriangular[Range@Binomial[1000, 2], 1000]; // Timing


{9.875, Null}

## Edit

Using the PadLeft ,come from @Algohi's answer

upperTriangular2[expr_?VectorQ, n_Integer] /; Length@expr == (n (n - 1))/2 :=
PadLeft[#, n] & /@ (Append[
Complement[#2, #1] & @@@
Partition[
Prepend[Take[expr, #] & /@ Accumulate@Reverse@Range[n - 1], {}],2, 1], {}])


Using TakeList:

Clear["Global*"];
n = 8;
list = Range[Binomial[n, 2]];

TakeList[list, Range[0, n - 1]] // Map[PadRight[#, n] &] // MatrixForm
TakeList[list, Range[n - 1, 0, -1]] //
Map[PadLeft[#, n] &] // MatrixForm


Using ReplacePart:

Clear["Global*"];
n = 8;
list = Range[Binomial[n, 2]];
posu = Flatten[#, 1] &@Table[{i, j}, {i, 1, n}, {j, i + 1, n}];
rulesup = Thread[posu -> list];
ReplacePart[ConstantArray[0, {n, n}], rulesup] // MatrixForm

posdn = Flatten[#, 1] &@Table[{i, j}, {i, 2, n}, {j, 1, i - 1}];
rulesdn = Thread[posdn -> list];
ReplacePart[ConstantArray[0, {n, n}], rulesdn] // MatrixForm
`

(* same outputs *)