Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
    
@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. –  Szabolcs Jul 23 at 14:04
    
The elements are machine integers only? Or...? –  Oleksandr R. Jul 23 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. –  Szabolcs Jul 23 at 14:16
    
Sorry, I just read "not to bias the answers" -- would you prefer that I delete the copied answer I just posted? –  Mr.Wizard Jul 23 at 14:59
    
@Mr.Wizard No, keep it! –  Szabolcs Jul 23 at 17:29

7 Answers 7

up vote 6 down vote accepted

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:

Developer`PackedArrayQ @ %

(* 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:

CompiledFunctionTools`CompilePrint@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
]
share|improve this answer
    
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." –  Tangshutao Sep 13 at 9:32
    
@Tangshutao I have added a note concerning this message to the bottom of my response. –  WReach Sep 13 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] –  Tangshutao Sep 14 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. –  WReach Sep 14 at 13:31

Probably not fast, but very simple:

upperTriangular[elements_, n_] :=
 SparseArray[
  Thread[
    SymmetrizedIndependentComponents[{n, n}, Antisymmetric[{1, 2}]] ->    
    elements]]
share|improve this answer
2  
SymmetrizedIndependentComponents? These function names are getting out of hand. +1 nevertheless as I don't recall seeing that function used before. –  Mr.Wizard Jul 23 at 17:11
    
Great! I didn't mention it but I was hoping for something using this functionality. Is StructuredArray useful here? –  Szabolcs Jul 23 at 17:25

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]
share|improve this answer

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} *)
share|improve this answer
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
share|improve this answer
    
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] –  Tangshutao Sep 14 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], {}])
share|improve this answer
ClearAll[f1];
f1 = RotateLeft[PadLeft[Internal`PartitionRagged[Range[# (#-1)/2],Range[#-1, 1, -1]], {#, #}]] &;

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

enter image description here

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

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

enter image description here

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

enter image description here

share|improve this answer
    
Your solution is very efficient –  Tangshutao Sep 13 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:) –  kguler Sep 13 at 15:13
    
Dear kguler, How can I understand Internal-PartitionRagged, I cannot find the usage of this function in Wolfram Help Documentation –  Tangshutao Sep 14 at 6:21
    
@Tangshuato, ??Internal`PartitionRagged gives some info about the function. To see the code you can use ClearAttributes[Internal`PartitionRagged, ReadProtected] followed by ??Internal`PartitionRagged. –  kguler Sep 14 at 6:27

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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