Doing symbolic algebra with matrix multiplication is very tedious, so I tried to do that in Julia. I take the following example concerning changing variance components from sire-mgs to calf-maternal (Kriese et al., 1991 and other places):
I use Symbolics.jl as follows:
> @variables s2s ssm s2m
G0=[s2s ssm; ssm s2m]
> L*G0*L'
> simplify.(L*G0*L',expand=true)
2×2 Matrix{Num}:
4s2s -2s2s + 4ssm
-2s2s + 4ssm 4s2m + s2s - 4ssm
and I get same as above. Now I want to get in the more classical form where variance components are stack into a vector (typically the upper or lower triangular of G0). We can use some "vec" algebra (from Henderson and Searle, http://www.jstor.org/stable/3315017 ) to do:
> kron(L,L)*vec(G0)
4-element Vector{Num}:
4s2s
-2s2s + 4ssm
-2s2s + 4ssm
4s2m + s2s - 4ssm
> vec(G0)
4-element Vector{Num}:
s2s
ssm
ssm
s2m
> kron(L,L)
4×4 Matrix{Int64}:
4 0 0 0
-2 4 0 0
-2 0 4 0
1 -2 -2 4
Coda: if we want to use vech, we have this (from google AI):
function vechh(A::AbstractMatrix{T}) where T
# Ensure matrix is square
m, n = size(A)
m == n || throw(DimensionMismatch("Matrix must be square"))
# Pre-allocate output vector of size n*(n+1)/2
len = (n * (n + 1)) >> 1
v = Vector{T}(undef, len)
k = 0
# Column-major order loops (j first, then i)
for j in 1:n
for i in j:m
@inbounds v[k += 1] = A[i, j]
end
end
return v
end
