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 $Aotimes Bto Botimes 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)

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

vech[W_] := With[{n = Length[W]}, 
  Flatten[MapThread[#1[[-#2 ;;]] &, {Transpose[W], Reverse@Range[n]}]]]

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 *)

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

vec[A]
vec[(A + Transpose[A])/2]
vec[As]