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;

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}];

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)

enter image description here

enter image description here

  • 1
    $\begingroup$ perhaps kMat[n_, m_] := IdentityMatrix[n m ][[Flatten[Range[Range[n], n m, n]]]] for Kmat? $\endgroup$
    – kglr
    Sep 15 '20 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]}, 
  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 *)
  • $\begingroup$ For $N_n$ it gets slow for $d \ge 5$ unfortunately. $\endgroup$
    – flinty
    Sep 15 '20 at 21:54
  • $\begingroup$ ^ 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] $\endgroup$
    – flinty
    Sep 15 '20 at 22:07
  • $\begingroup$ 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 $\endgroup$ Sep 15 '20 at 22:11
  • $\begingroup$ There's definitely a pattern to $G_n$ which can be seen in gdupe[6] // 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. $\endgroup$
    – flinty
    Sep 15 '20 at 22:13
  • 1
    $\begingroup$ @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. $\endgroup$
    – flinty
    Sep 15 '20 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;

with this your examples read:

vec[(A + Transpose[A])/2]
  • $\begingroup$ You need vech operator -- note that duplication matrix G_n is not square $\endgroup$ Sep 15 '20 at 20:08
  • $\begingroup$ vech[] may be coded as: vech[m_] := Flatten@Table[m[[i, j]], {i, Length[m]}, {j, i, Length[m[[1]]]}] $\endgroup$ Sep 16 '20 at 8:22

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.