4
$\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
    $\begingroup$ perhaps kMat[n_, m_] := IdentityMatrix[n m ][[Flatten[Range[Range[n], n m, n]]]] for Kmat? $\endgroup$ – kglr Sep 15 at 18:32
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 *)
| improve this answer | |
$\endgroup$
  • $\begingroup$ For $N_n$ it gets slow for $d \ge 5$ unfortunately. $\endgroup$ – flinty Sep 15 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 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$ – Yaroslav Bulatov Sep 15 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 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 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] 
| improve this answer | |
$\endgroup$
  • $\begingroup$ You need vech operator -- note that duplication matrix G_n is not square $\endgroup$ – Yaroslav Bulatov Sep 15 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$ – Daniel Huber Sep 16 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.