Sorry for my late answer: I had several deadlines and I couldn’t work on this item.
Luckily, now I had a little time.
Mathematical Background
First of all, I want to link some resources
1 The elimination Matrix
2 Symmetry, 0-1 Matrices and Jacobians
Given 4 matrices A, B, C, D, the following identity holds
(\operatorname{vec} A)^{\prime}(E \otimes B) \operatorname{vec} C=\operatorname{tr} A^{\prime} BCE^{\prime}\qquad (1)
if the right-hand side exists.
If A is a m,n matrix, then \operatorname{vec} A is defined by
\operatorname{vec} A=\left(\begin{array}{c}
A_{\cdot 1} \\
\vdots \\
A_{\cdot n}
\end{array}\right)
If A is a square symmetric matrix, \operatorname{vecp} A is the collection of all unique elements of A.
\operatorname{vec} A and \operatorname{vecp} A can be mapped one to another using the duplication D and elimination L (their definitions can be found in references 1 and 2; in particular, reference 1 gives an explicit definition of the elimination and duplication matrix):
\operatorname{vec} A = D \operatorname{vecp} A \qquad (2)
\operatorname{vecp} A = L\operatorname{vec} A
In particular, Equation (1) can be rewritten using Equation (2)
(\operatorname{vecp} A)^{\prime} D^{\prime}(E \otimes B) D\operatorname{vecp} C=\operatorname{tr} A^{\prime} BCE^{\prime}\qquad (3)
Implementation
Here is my implementation of the duplication and elimination matrix
function uᵢⱼ(n::Int, i::Int, j::Int)
u = zeros(floor(Int,0.5*n*(n+1)))
u[floor(Int64, (j-1)*n+i-0.5*j*(j-1))] = 1
return u
end
function eᵢ(n::Int, i::Int)
e = zeros(n)
e[i] = 1
return e
end
function Eᵢⱼ(n::Int, i::Int, j::Int)
E = zeros((n,n))
E[i,j] = 1
return E
end
function Tᵢⱼ(n::Int, i::Int, j::Int)
if i != j
T = Eᵢⱼ(n,i,j) + Eᵢⱼ(n,j,i)
else
T = Eᵢⱼ(n,i,i)
end
return T
end
function EliminationMatrix(n::Int)
L = zeros((floor(Int64,0.5*n*(n+1)), n*n))
for i in 1:n
for j in 1:i
L += kron(uᵢⱼ(n,i,j), transpose(vec(Eᵢⱼ(n,i,j))))
end
end
return L
end
function DuplicationMatrix(n::Int)
D = zeros((n*n, floor(Int64,0.5*n*(n+1))))
for i in 1:n
for j in 1:i
D += transpose(kron(uᵢⱼ(n,i,j), transpose(vec(Tᵢⱼ(n,i,j)))))
end
end
return D
end
Checks
Using some random matrix, I checked that the algebric properties are satisfied
#generating random matrices
n = 10
A = rand(n,n)
A = A + transpose(A)
B = rand(n,n)
B = B + transpose(B)
C = rand(n,n)
C = C + transpose(C)
E = rand(n,n)
E = E + transpose(E)
AB = zeros(n,n)
CE = zeros(n,n)
ABCE = zeros(n,n)
LinearAlgebra.mul!(AB, transpose(A), B)
LinearAlgebra.mul!(CE, C, transpose(E))
LinearAlgebra.mul!(ABCE, AB, CE)
#result with first technique
LinearAlgebra.tr(ABCE)
kronEB = kron(E,B)
Dáµ€kronEB = zeros(floor(Int,n*0.5*(n+1)),n*n)
LinearAlgebra.mul!(Dáµ€kronEB, Transpose(DuplicationMatrix(n)), kronEB)
Dáµ€kronEBD = zeros(floor(Int,n*0.5*(n+1)),floor(Int,n*0.5*(n+1)))
LinearAlgebra.mul!(Dáµ€kronEBD, Dáµ€kronEB, (DuplicationMatrix(n)))
vecpA = zeros(floor(Int,n*0.5*(n+1)))
LinearAlgebra.mul!(vecpA, EliminationMatrix(n), vec(A))
vecpC = zeros(floor(Int,n*0.5*(n+1)))
LinearAlgebra.mul!(vecpC, EliminationMatrix(n), vec(C))
vecpAáµ€Dáµ€kronEBD = zeros(1,floor(Int,n*0.5*(n+1)))
LinearAlgebra.mul!(vecpAáµ€Dáµ€kronEBD, transpose(vecpA), Dáµ€kronEBD)
vecpAáµ€Dáµ€kronEBDvecpC = zeros(1,1)
#result with the second method
LinearAlgebra.mul!(vecpAáµ€Dáµ€kronEBDvecpC , vecpAáµ€Dáµ€kronEBD, vecpC)
The results agree (up to machine precision).
@mcabbott , are you aware of any packages including the duplication and elimination matrix as implemented here? I looked into LinearAlgebra, but I found nothing similar.
Thank you,
Marco