Duplication and elimination matrices

From HandWiki

In mathematics, especially in linear algebra and matrix theory, the duplication matrix and the elimination matrix are linear transformations used for transforming half-vectorizations of matrices into vectorizations or (respectively) vice versa.

Duplication matrix

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

[math]\displaystyle{ D_n \mathrm{vech}(A) = \mathrm{vec}(A) }[/math].

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

[math]\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} }[/math]


The explicit formula for calculating the duplication matrix for a [math]\displaystyle{ n \times n }[/math] matrix is:

[math]\displaystyle{ D^T_n = \sum \limits_{i \ge j} u_{ij} (\mathrm{vec}T_{ij})^T }[/math]

Where:

  • [math]\displaystyle{ u_{ij} }[/math] is a unit vector of order [math]\displaystyle{ \frac{1}{2} n (n+1) }[/math] having the value [math]\displaystyle{ 1 }[/math] in the position [math]\displaystyle{ (j-1)n+i - \frac{1}{2}j(j-1) }[/math] and 0 elsewhere;
  • [math]\displaystyle{ T_{ij} }[/math] is a [math]\displaystyle{ n \times n }[/math] matrix with 1 in position [math]\displaystyle{ (i,j) }[/math] and [math]\displaystyle{ (j,i) }[/math] 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();
}

Elimination matrix

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

[math]\displaystyle{ L_n \mathrm{vec}(A) = \mathrm{vech}(A) }[/math][1]

By the explicit (constructive) definition given by (Magnus Neudecker), the [math]\displaystyle{ \frac{1}{2}n(n+1) }[/math] by [math]\displaystyle{ n^2 }[/math] elimination matrix [math]\displaystyle{ L_n }[/math] is given by

[math]\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), }[/math]

where [math]\displaystyle{ e_i }[/math] is a unit vector whose [math]\displaystyle{ i }[/math]-th element is one and zeros elsewhere, and [math]\displaystyle{ E_{ij} = e_ie_j^T }[/math].

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 [math]\displaystyle{ 2 \times 2 }[/math] matrix [math]\displaystyle{ A = \left[\begin{smallmatrix} a & b \\ c & d \end{smallmatrix}\right] }[/math], one choice for this transformation is given by

[math]\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} }[/math].

Notes

  1. (Magnus Neudecker), Definition 3.1

References

de:Eliminationsmatrix