symbolic algebra in Julia

 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):

(this is from an unpublished note from Peter Sullivan [Lactanet] and Urs Schuler [Qualitas])

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:

which we do as:

> kron(L,L)*vec(G0)

4-element Vector{Num}:

              4s2s

      -2s2s + 4ssm

      -2s2s + 4ssm

 4s2m + s2s - 4ssm

or, if we only want the coefficients:

> 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