13
$\begingroup$

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.

Improve this question
$\endgroup$
3
  • $\begingroup$ @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. $\endgroup$
    – Szabolcs
    Jul 23, 2014 at 14:04
  • 1
    $\begingroup$ The elements are machine integers only? Or...? $\endgroup$
    – Oleksandr R.
    Jul 23, 2014 at 14:14
  • $\begingroup$ @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. $\endgroup$
    – Szabolcs
    Jul 23, 2014 at 14:16

7 Answers 7

Reset to default
10
$\begingroup$

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
]
Improve this answer
$\endgroup$
4
  • $\begingroup$ 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." $\endgroup$
    – xyz
    Sep 13, 2014 at 9:32
  • $\begingroup$ @Tangshutao I have added a note concerning this message to the bottom of my response. $\endgroup$
    – WReach
    Sep 13, 2014 at 12:07
  • $\begingroup$ 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] $\endgroup$
    – xyz
    Sep 14, 2014 at 6:41
  • $\begingroup$ @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. $\endgroup$
    – WReach
    Sep 14, 2014 at 13:31
12
$\begingroup$

Probably not fast, but very simple:

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

Update: Recently bumped into Statistics`Library`VectorToUpperTriangularMatrix 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 = Statistics`Library`VectorToUpperTriangularMatrix[#, #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)$

Original 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

Improve this answer
$\endgroup$
4
  • $\begingroup$ Your solution is very efficient $\endgroup$
    – xyz
    Sep 13, 2014 at 14:47
  • $\begingroup$ @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:) $\endgroup$
    – kglr
    Sep 13, 2014 at 15:13
  • $\begingroup$ Dear kguler, How can I understand Internal-PartitionRagged, I cannot find the usage of this function in Wolfram Help Documentation $\endgroup$
    – xyz
    Sep 14, 2014 at 6:21
  • $\begingroup$ @Tangshuato, ??InternalPartitionRagged_ gives some info about the function. To see the code you can use _ClearAttributes[InternalPartitionRagged, ReadProtected] followed by ??Internal`PartitionRagged. $\endgroup$
    – kglr
    Sep 14, 2014 at 6:27
6
$\begingroup$

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} *)
Improve this answer
$\endgroup$
6
$\begingroup$

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]
Improve this answer
$\endgroup$
4
$\begingroup$ 
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
Improve this answer
$\endgroup$
1
  • $\begingroup$ 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] $\endgroup$
    – xyz
    Sep 14, 2014 at 6:01
1
$\begingroup$

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], {}])
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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