# Commutation, symmetrizer and duplication matrices

The following 3 matrices are useful when viewing matrices as vectors, known as commutation $$K_n$$, symmetrizer $$N_n$$ and duplication $$G_n$$. They are usually defined by their matrix relations below.

$$\begin{eqnarray} \text{vec}A & = & K_n \text{vec}A' \\ \text{vec}((A+A')/2) & = &N_n \text{vec}A\\ \text{vec}A_s & = & G_n \text{vech}A_s\\ \end{eqnarray}$$

Here $$\text{vec}$$ is a vectorization operator that stacks columns, and $$\text{vech}$$ is "lower-half" vectorization, stacking columns of the lower half of the matrix. $$A$$ is arbitrary matrix, $$A_s$$ is symmetric

(A related matrix commutes the order of Kronecker product $$A\otimes B\to B\otimes A$$)

I have an ugly-looking implementation of the first two matrices based on some some algebra done by Seber, "Handbook of Statistics", section 11.5. Can someone see a good way to implement the third matrix?

Also wondering if there's some functionality in Mathematica that would obviate the need to do manual algebra and instead rely on matrix relations above.

(* Commutation matrix m,n *)

Kmat[m_, n_] := Module[{x, X, before, after, positions, matrix},
X = Array[x, {m, n}];
before = Flatten@vec@X;
after = Flatten@vec@Transpose[X];
positions =
MapIndexed[{First@#2, First@Flatten@Position[before, #]} &, after];
matrix = SparseArray[# -> 1 & /@ positions] // Normal
];

Nmat[n_] := (Normal@Kmat[n, n] + IdentityMatrix[n^2])/2;
Gmat[n_] := Array[1 &, {n, n (n + 1)/2}];

n = 3;

Clear[a];
A = Array[a, {3, 3}];
As = Array[a[Min[#1, #2], Max[#1, #2]] &, {n, n}];

vec[W_] := Transpose@{Flatten@Transpose[W]};
vech[W_] := Flatten@Table[Table[W[[i, j]], {i, j, n }], {j, 1, n}];

On[Assert];
Assert[vec[A] == Kmat[n, n].vec[A\[Transpose]]]
Assert[vec[(A + Transpose[A])/2] == Nmat[n].vec[A] // Reduce]
Assert[vec[As] == Gmat[n].vech[As] // Reduce]


# Official description

Here's description from Seber's Handbook of Statistics: ($$G_3=D_3$$ is duplication matrix, $$H_3$$ is it's inverse -- the elimination matrix, and $$I_{(3,3)}$$ is the commutation matrix)  • perhaps kMat[n_, m_] := IdentityMatrix[n m ][[Flatten[Range[Range[n], n m, n]]]] for Kmat? – kglr Sep 15 at 18:32

I hope this does the trick. It's more code than yours but I've come at it from a slightly different angle - I suppose another implementation can't hurt right? I've used FindPermutation to get $$K_n$$ and SolveAlways for non-square $$G_n$$:

vec[W_] := Join @@ Transpose[W]

vech[W_] := With[{n = Length[W]},

getperm[perm_, n_] := Permute[IdentityMatrix[n*n], perm]

kcomm[n_] := With[{mtx = ArrayReshape[Range[n*n], {n, n}]},
getperm[FindPermutation[vec[Transpose[mtx]], vec[mtx]], Length[mtx]]]

nsymm[n_] := (kcomm[n] + IdentityMatrix[n^2])/2

gdupe[n_] :=
With[{mtx = Array[a[Min[#1, #2], Max[#1, #2]] &, {n, n}],
gmatrix = Array[x, {n*n, n (n + 1)/2}]},
gmatrix /. First[SolveAlways[vec[mtx] == gmatrix.vech[mtx], Variables[mtx]]]]

(* tests *)
d = 3;
m = RandomReal[{-1, 1}, {d, d}];
kcomm[d].vec[Transpose[m]] == vec[m]
(* True *)

nsymm[d].vec[m] == vec[(m + Transpose[m])/2]
(* True *)

vec[Normal[Symmetrize[m]]] == gdupe[d].vech[Normal[Symmetrize[m]]]
(* True *)

• For $N_n$ it gets slow for $d \ge 5$ unfortunately. – flinty Sep 15 at 21:54
• ^ just saw your version which uses $K_n$ to generate $N_n$ so I guess we can do this which is much faster: With[{n = 5}, (Normal@kcomm[n] + IdentityMatrix[n^2])/2] – flinty Sep 15 at 22:07
• Nice trick with SolveAlways! It's expected that it won't be as fast as the manually coded algebraic solution, but good trick to have in one's toolbox – Yaroslav Bulatov Sep 15 at 22:11
• There's definitely a pattern to $G_n$ which can be seen in gdupe // MatrixPlot. 6 on the diagonal, then five but shifted, then four, etc.. with various interruptions. Perhaps there's another way to construct it which is much faster. – flinty Sep 15 at 22:13
• @YaroslavBulatov there is apparently a form here which involves a sum en.wikipedia.org/wiki/Duplication_and_elimination_matrices , but I can't quite get it to work. – flinty Sep 15 at 23:01

If I understand correctly, then you only need the operator "vec". This is clear for the first line. The second line applies vec to the symmetrized version of A: (A+Transpose[A])/2. And the third line applies "vec" to a symmetric matrix, the operator is the same, only the operand is different. Therefor in MMA I would code:

A = Array[a, {3, 3}];
As = Array[a[Min[#1, #2], Max[#1, #2]] &, {n, n}];
vec[m_]:= List /@ Flatten@Transpose@m;


vec[A]