The duplication matrix ${\displaystyle D_{n}}$  is the unique ${\displaystyle n^{2}\times {\frac {n(n+1)}{2}}}$  matrix which, for any ${\displaystyle n\times n}$  symmetric matrix ${\displaystyle A}$ , transforms ${\displaystyle \mathrm {vech} (A)}$  into ${\displaystyle \mathrm {vec} (A)}$ :

${\displaystyle D_{n}\mathrm {vech} (A)=\mathrm {vec} (A)}$ .

For the ${\displaystyle 2\times 2}$  symmetric matrix ${\displaystyle A=\left[{\begin{smallmatrix}a&b\\b&d\end{smallmatrix}}\right]}$ , this transformation reads

${\displaystyle D_{n}\mathrm {vech} (A)=\mathrm {vec} (A)\implies {\begin{bmatrix}1&0&0\\0&1&0\\0&1&0\\0&0&1\end{bmatrix}}{\begin{bmatrix}a\\b\\d\end{bmatrix}}={\begin{bmatrix}a\\b\\b\\d\end{bmatrix}}}$ 

The explicit formula for calculating the duplication matrix for a ${\displaystyle n\times n}$  matrix is:

${\displaystyle D_{n}^{T}=\sum \limits _{i\geq j}u_{ij}(\mathrm {vec} T_{ij})^{T}}$ 

Where:

• ${\displaystyle u_{ij}}$  is a unit vector of order ${\displaystyle {\frac {1}{2}}n(n+1)}$  having the value ${\displaystyle 1}$  in the position ${\displaystyle (j-1)n+i-{\frac {1}{2}}j(j-1)}$  and 0 elsewhere;
• ${\displaystyle T_{ij}}$  is a ${\displaystyle n\times n}$  matrix with 1 in position ${\displaystyle (i,j)}$  and ${\displaystyle (j,i)}$  and zero elsewhere

Here is a C++ function using Armadillo (C++ library):

arma::mat duplication_matrix(const int &n) {
    arma::mat out((n*(n+1))/2, n*n, arma::fill::zeros);
    for (int j = 0; j < n; ++j) {
        for (int i = j; i < n; ++i) {
            arma::vec u((n*(n+1))/2, arma::fill::zeros);
            u(j*n+i-((j+1)*j)/2) = 1.0;
            arma::mat T(n,n, arma::fill::zeros);
            T(i,j) = 1.0;
            T(j,i) = 1.0;
            out += u * arma::trans(arma::vectorise(T));
        }
    }
    return out.t();
}

An elimination matrix ${\displaystyle L_{n}}$  is a ${\displaystyle {\frac {n(n+1)}{2}}\times n^{2}}$  matrix which, for any ${\displaystyle n\times n}$  matrix ${\displaystyle A}$ , transforms ${\displaystyle \mathrm {vec} (A)}$  into ${\displaystyle \mathrm {vech} (A)}$ :

${\displaystyle L_{n}\mathrm {vec} (A)=\mathrm {vech} (A)}$ . ^[1]

By the explicit (constructive) definition given by Magnus & Neudecker (1980), the ${\displaystyle {\frac {1}{2}}n(n+1)}$  by ${\displaystyle n^{2}}$  elimination matrix ${\displaystyle L_{n}}$  is given by

${\displaystyle L_{n}=\sum _{i\geq j}u_{ij}\mathrm {vec} (E_{ij})^{T}=\sum _{i\geq j}(u_{ij}\otimes e_{j}^{T}\otimes e_{i}^{T}),}$ 

where ${\displaystyle e_{i}}$  is a unit vector whose ${\displaystyle i}$ -th element is one and zeros elsewhere, and ${\displaystyle E_{ij}=e_{i}e_{j}^{T}}$ .

Here is a C++ function using Armadillo (C++ library):

arma::mat elimination_matrix(const int &n) {
    arma::mat out((n*(n+1))/2, n*n, arma::fill::zeros);
    for (int j = 0; j < n; ++j) {
        arma::rowvec e_j(n, arma::fill::zeros);
        e_j(j) = 1.0;
        for (int i = j; i < n; ++i) {
            arma::vec u((n*(n+1))/2, arma::fill::zeros);
            u(j*n+i-((j+1)*j)/2) = 1.0;
            arma::rowvec e_i(n, arma::fill::zeros);
            e_i(i) = 1.0;
            out += arma::kron(u, arma::kron(e_j, e_i));
        }
    }
    return out;
}

For the ${\displaystyle 2\times 2}$  matrix ${\displaystyle A=\left[{\begin{smallmatrix}a&b\\c&d\end{smallmatrix}}\right]}$ , one choice for this transformation is given by

${\displaystyle L_{n}\mathrm {vec} (A)=\mathrm {vech} (A)\implies {\begin{bmatrix}1&0&0&0\\0&1&0&0\\0&0&0&1\end{bmatrix}}{\begin{bmatrix}a\\c\\b\\d\end{bmatrix}}={\begin{bmatrix}a\\c\\d\end{bmatrix}}}$ .

• Magnus, Jan R.; Neudecker, Heinz (1980), "The elimination matrix: some lemmas and applications", SIAM Journal on Algebraic and Discrete Methods, 1 (4): 422–449, doi:10.1137/0601049, ISSN 0196-5212.
• Jan R. Magnus and Heinz Neudecker (1988), Matrix Differential Calculus with Applications in Statistics and Econometrics, Wiley. ISBN 0-471-98633-X.
• Jan R. Magnus (1988), Linear Structures, Oxford University Press. ISBN 0-19-520655-X