7
$\begingroup$

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)

enter image description here

enter image description here

$\endgroup$
1
  • 1
    $\begingroup$ perhaps kMat[n_, m_] := IdentityMatrix[n m ][[Flatten[Range[Range[n], n m, n]]]] for Kmat? $\endgroup$
    – kglr
    Sep 15, 2020 at 18:32

2 Answers 2

4
$\begingroup$

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 *)
$\endgroup$
6
  • $\begingroup$ For $N_n$ it gets slow for $d \ge 5$ unfortunately. $\endgroup$
    – flinty
    Sep 15, 2020 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, 2020 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, 2020 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, 2020 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, 2020 at 23:01
0
$\begingroup$

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]
vec[(A + Transpose[A])/2]
vec[As] 
$\endgroup$
2
  • $\begingroup$ You need vech operator -- note that duplication matrix G_n is not square $\endgroup$ Sep 15, 2020 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, 2020 at 8:22

This site is temporarily in read-only mode and not accepting new answers.

Not the answer you're looking for? Browse other questions tagged .