# Efficient parametrization of the covariance matrix with some covariances constrained to zero

I'm trying to estimate the unknown 8x8 covariance matrix X in R using the maximum likehood, but I have problems of figuring out the efficient way of parametrization of X when some of the covariances in X are constrained to zero. When there's no constraints, I have used a cholesky factorization, ie. I have parameters which correspond to the lower triangular matrix L, and I get X by X<-t(L)%*%L which is a proper positive definite matrix, but how to do it now that I want some of the cells of X to be constrained to zero? In my case, X is constrained in a way that the upper left and lower right 4x4 matrices can be anything (but of course having the properties of proper covariance matrices), and the off-diagonal 4x4 matrix is a diagonal matrix, ie. the the first variable x1 correlates with variables x2,x3,x4 and x5, variable x2 correlates with x1,x3,x4 and x6 etc.

Thanks.

You could start from your original approach and impose the equations that the specified coefficients are 0. This leads to a fairly large system of polynomial equations - 12 equations (corresponding to the zeroes in the upper triangle) in 25 unknowns. I can solve this in Maple and obtain a union of 87 parametric solutions of dimensions from 11 to 16, each of which parametrizes a subset of the full solution and the union of which forms that full solution.

For example, one solution (one of the two 16-dimensional ones) is given by:

$$\begin{gather} A_{{1,1}}=0,A_{{2,1}}=0,A_{{2,2}}=0,A_{{3,1}}=0,A_{{3,2}}=0,A _{{3,3}}=0, \\ A_{{5,1}}=-{\frac {A_{{4,2}}A_{{5,2}}+A_{{4,3}}A_{{5,3}}+A_ {{4,4}}A_{{5,4}}}{A_{{4,1}}}}, \\ A_{{6,1}}=-{\frac {A_{{4,2}}A_{{6,2}}+A_ {{4,3}}A_{{6,3}}+A_{{4,4}}A_{{6,4}}}{A_{{4,1}}}}, \\ A_{{7,1}}=-{\frac {A_ {{4,2}}A_{{7,2}}+A_{{4,3}}A_{{7,3}}+A_{{4,4}}A_{{7,4}}}{A_{{4,1}}}}. \end{gather}$$

However, even though all of these solutions are needed to describe all lower-triangular matrices for which $L L^T$ is of the form you require, it may well be that to just obtain all such matrices, we can make do with fewer solutions. In fact, I have a sneaking suspicion that the solution above might cover all cases. This is just a hunch that would require a more thorough investigation, though.

In order to reproduce the computation if you have a copy of Maple (15, in my case) available, you can run the following:

# Construct a generic lower triangular 8x8 matrix.
lt := Matrix(8, symbol = A, shape = triangular[lower]):
mm := lt . lt^%T:

# Select the positions that should be zero.
positions := [seq(seq([i, j], j = 5 .. 8), i = 1 .. 4)]:
positions := remove(pair -> pair[2] = pair[1] + 4, positions):

# Construct and solve the system of equations.
sys := map(pair -> mm[op(pair)], positions):
solutions := [solve(sys)]:
nops(solutions); # returns 87, so 87 different solutions.

# Split each solution into trivial equations (like A[2,3] = A[2,3]) and
# nontrivial ones. This separates free variables from those determined
# by other ones.
split := map2([selectremove], evalb, solutions):
numbers := map2(map, nops, split):
convert(numbers, set); # returns {[11, 14], [12, 13], ..., [16, 9]},
# showing that the dimension of the components runs
# from 11 to 16.

# Select the solutions of dimension 16.
dim16 := select(pair -> nops(pair[1]) = 16, split):
dim16 := map(pair -> pair[2], dim16): # We're only interested in the
# nontrivial equations.
nops(dim16); # returns 2, showing there are 2. One of these is the one
# printed above.