# R package to solve Gaussian MLE under conditional independence constraints

Is there any R package or function to solve Gaussian MLE under conditional independence constraints?

Suppose we have $$y_i\overset{i.i.d}{\sim}\mathcal{N}(0,\Sigma_{p\times p})$$, $$i = 1,2,\ldots,n$$. We know that $$(\Sigma^{-1})_{ij} = 0$$, for some $$i,j\in\{1,2,\ldots,p\}$$, i.e. $$X_i$$ and $$X_j$$ are conditionally independent given the rest of the variables. We would like to find MLE of $$\Sigma$$ under the conditional independence constraints.

I also tried to implement it according to section 5.1 of the above paper, but I couldn't do it successfully.

I wonder if there is any R implementation to find Gaussian MLE under conditional independence constraints?

I finished the implementation myself. Here is the code.

gaussian_MLE = function(S, omega_structure, initial_omega = initial_omega, maxit = 2e9, tol = 1e-15){
omegas = list()

# S is sample covariance matrix; omega_structure is the adjacency matrix.

temp = which(omega_structure == 1, arr.ind = TRUE)
temp = temp[temp[,1]>=temp[,2],]

r = 1
converged = FALSE

omega_old = initial_omega
omega_old[!omega_structure] = 0
omega_current = omega_old
while(!converged){

omega_old = omega_current

for(i in 1:nrow(temp)){

S_current = solve(omega_current)

u = temp[i,1]
v = temp[i,2]

sqrt = sqrt(S_current[u,u] * S_current[v,v])

rho_bar = S[u, v]/sqrt
# print(rho_bar)

rho = S_current[u, v]/sqrt

if(rho==1){
if(rho_bar > 0){
s = (1-rho_bar)/(2*rho_bar)/sqrt
}else{
stop("here1")
}
}else if(rho == -1){
if(rho_bar < 0){
s = (1+rho_bar)/(2*rho_bar)/sqrt
}else{
stop("here2")
}
}else{
one_minus_rho_square = 1-rho^2
a = rho_bar*one_minus_rho_square
b = -(one_minus_rho_square+2*rho*rho_bar)
c = rho-rho_bar

if(rho_bar==0){
s = rho/one_minus_rho_square/sqrt
}else if(rho_bar>0){
s = min((-b+sqrt(b^2-4*a*c))/(2*a)/sqrt,(-b-sqrt(b^2-4*a*c))/(2*a)/sqrt)
}else{
s = max((-b+sqrt(b^2-4*a*c))/(2*a)/sqrt,(-b-sqrt(b^2-4*a*c))/(2*a)/sqrt)
}

}

if(u==v){
omega_current[u,v] =  omega_current[u,v] +  s
}else{
omega_current[u,v] =  omega_current[u,v] +  s
omega_current[v,u] =  omega_current[v,u] +  s
}

}

if(max(abs(omega_current - omega_old))<tol){
converged = TRUE
}else{
r = r+1
if((r %% 1000000)==0){
print(r)
}
}

if(r>maxit){
stop("Algorithm did not converge.")
break;
}

# print(r)
# print(omega_old)

# omegas[[length(omegas)+1]] = omega_old
# omegas <<- omegas

}

if(max(abs(solve(omega_current)*omega_structure - S*omega_structure))>1e-13){
stop("here3")
}
eigen(omega_current)

rrr <<- r

return(omega_current)

}