How to implement a Gaussian NN with covariance matrix?

That’s probably the way to go. Just wanted to chime in to correct a minor mistake in c2idx:
Looking at the vector indices this translates to:

julia> c2idx(i, j) = (i*(i-1))÷2+j
c2idx (generic function with 1 method)

julia> [ifelse(j<i, 0, c2idx(i,j)) for i in 1:N, j in 1:N]
4×4 Matrix{Int64}:
 1  2  3   4
 0  3  4   5
 0  0  6   7
 0  0  0  10

i think you missed the sidelength in the formula:

julia> new_c2idx(i, j) = ((2N-i)*(i-1))÷2+j
new_c2idx (generic function with 1 method)

julia> [ifelse(j<i, 0, new_c2idx(i,j)) for i in 1:N, j in 1:N]
4×4 Matrix{Int64}:
 1  2  3   4
 0  5  6   7
 0  0  8   9
 0  0  0  10
2 Likes