Python code for Duplication matrices

Question

I was looking for a resource/code/script in python for forming duplication matrices but surprisingly it still evades me. The wiki page gives a pseudo-code for generating this 'unique' duplication matrix of size $N^2 \times N(N+1)/2$ for a symmetric matrix of size $N$.

Basically, I want matrix $D$ such that for an symmetric matrix $A \in \mathbb{R}^{N \times N}$,

$$D \: \text{vech}(A) = \text{vec}(A)$$

where vec and vech denote vectorization and half-vectorization, respectively.

Here is what I could manage from the wiki page

import numpy as np
def duplication_matrix(N):
    # source: https://en.wikipedia.org/wiki/Duplication_and_elimination_matrices
    N_bar = int(N * (N + 1) / 2.)
    # duplication matrix initialization
    D = np.zeros((N ** 2, N_bar))

    for ii in range(N):
        for jj in range(N):
            if jj <= ii:
                u_ij = np.zeros((N_bar, 1))
                tmp_idx = int(jj * N + ii + 1. - ((jj + 1.) * jj) / 2.) - 1
                u_ij[tmp_idx] = 1.
                T_ij = np.zeros((N, N))
                T_ij[ii, jj] = 1.
                T_ij[jj, ii] = 1.

                tmp = u_ij.dot(vectorize(T_ij).reshape((N ** 2, 1)).transpose())
                D += tmp.transpose()

    return D

You might find the indexing a bit difference corresponding to variable tmp_idx compared to what is given on the wiki page but it is done just to account for the indexing in python which starts from $0$.

Problem is that it apparently worked for $N=2$. I tried it with $N = 4$ but it didn't yield the right result. Before I dive into creating something by myself, I wanted to know if someone has already done it as it is surprising to me that this is not readily available.

Edit: I was trying to do something proper from this post but there were inconsistencies with the notation. So, I have created a hackey answer that solves my problem. But I will still appreciate a proper answer!

Answer 1

This is a hackey way of doing this, but basically, I created a list of all alphabets (as I needed something unique to keep track of entries) and created a chararray of size $N \times N$. Then just by comparing the entries of the vectorized and half-vectorized entities, I can create the duplication matrix, $D$. Here it goes:

import numpy as np
def duplication_matrix_char(N):
    N_bar = int(N * (N + 1) / 2.)

    letters = ['a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n', 'o', 'p', 'q', 'r', 's', 't', 'u',
           'v', 'w', 'x', 'y', 'z']
    idx = np.triu_indices(N)

    # duplication matrix initialization
    D = np.zeros((N ** 2, N_bar))
    D_char = np.chararray((N, N))

    for ii in range(N_bar):
        D_char[idx[0][ii], idx[1][ii]] = letters[ii]
        D_char[idx[1][ii], idx[0][ii]] = letters[ii]

    # vectorizing D_char
    for kk in range(N):
        vecD_char[kk * N: (kk + 1) * N] = D_char[:, kk].reshape((N, 1))

    # half vectorize D_char
    vechD_char = np.chararray((N_bar, 1))
    idx_l = np.tril_indices(N)
    for mm in range(N_bar):
        vech[mm, 0] = D_char[idx_l[0][mm], idx_l[1][mm]]

    for jj in range(N ** 2):
        idx = np.where(vechD_char == vecD_char[jj])
        D[jj, idx[0]] = 1.

    return D

As there are $26$ alphabets, you can create till $N = 6$. You can easily extend this for any $N$ as long as you have $N(N+1)/2$ unique letters!

Answer 2

I worked on this years ago and came up with the following Julia function

function duplication(n)
   i = collect(1:n*n)
   J = zeros(Int, n,n)
   z = 1
   for x = 1:n
     for y = x:n
       J[x,y] = J[y,x] = z  # duplication occurs here
       z += 1
     end
   end
   # pass 3 vectors {i-index, j-index, value} to sparse()
   return sparse(i, vec(J), 1)
end

I also wrote a function based on the Wikipedia pseudo-code

function wiki_duplication(n)
   p = div(n*(n+1),2)
   A,D,u = zeros(Int,n,n), zeros(Int,n*n,p), zeros(Int,p,1)
   for i = 1:n
      for j = 1:i
         # use "." syntax to avoid reallocating the matrices
         A .= 0
         u .= 0
         ndx = i + (j-1)*n + div(j-j*j,2)
         A[i,j] = A[j,i] = u[ndx] = 1
         D += vec(A)*u'
      end
   end
   return D
end

This also works but is orders of magnitude slower than the first function.

Translating them into Python/NumPy shouldn't be terribly difficult.