Populate an upper triangular matrix from a vector of elements
Probably not fast, but very simple:
upperTriangular[elements_, n_] :=
SparseArray[
Thread[
SymmetrizedIndependentComponents[{n, n}, Antisymmetric[{1, 2}]] ->
elements]]
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
]
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]