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.
 A: 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.

