Estimating correlation matrix using numeric likelihood maximization I'm performing maximum likelihood estimation on jointly distributed data and I'm having some issues estimating the correlation terms. I am using an approach based on the Cholesky decomposition, but I feel that it has too many free parameters during estimation. Let me clarify the question.
Here's the formal definition of the data/model:
Suppose that $\mathbf{X}_i=(x_{i,1},x_{i,2},...,x_{i,n})$ and that $\mathbf{X}_i\sim Norm(\mathbf{0}_n,\mathbf{\Sigma})$, where $\mathbf{0}_n$ is just an $n$-long vector of zeros and $\mathbf{\Sigma}$ is the correlation (not covariance) matrix between the elements. Stated clearly:
$\mathbf{\Sigma}=\begin{bmatrix}
    1         & \rho_{12} & \rho_{13} & \dots  & \rho_{1n} \\
    \rho_{12} & 1         & \rho_{23} & \dots  & \rho_{2n} \\
    \vdots    & \vdots    & \vdots    & \ddots & \vdots \\
    \rho_{1n} & \rho_{2n} & \rho_{3n} & \dots  & 1
\end{bmatrix}$
Suppose further that all of my observations $\mathbf{X}_1,\mathbf{X}_2,...,\mathbf{X}_m$ are i.i.d.
What I'm trying to estimate are all the correlation terms contained in $\mathbf{\Sigma}$. Note that $\mathbf{\Sigma}$ contains $(n^2-n)/2$ correlation terms.
I know that I can't try to estimate the correlation terms directly, because during the estimation process that approach might generate correlation matrices that aren't positive semi-definite. What I am doing, however, is using the Cholesky decomposition approach.
Here, I generate a random vector of size $(n^2-n)/2+n$, place those elements in a lower triangular matrix, and multiply the result by its own transpose. Then, I normalize this new matrix by the using the square roots of the main diagonal. Formally:
$k = ((n^2-n)/2+n)$
$\mathbf{z} =(z_1, z_2, ..., z_k)=k $-$ long \ vector \ with \ iid \ random \ draws \ from \ any \ distribution$
$\mathbf{L} = \begin{bmatrix}
    z_1    & 0      & 0      & \dots  & 0 \\
    z_2    & z_3    & 0      & \dots  & 0 \\
    z_4    & z_5    & z_6    & \dots  & 0 \\
    \vdots & \vdots & \vdots & \ddots & \vdots \\
    z_{k-(n-1)} & z_{k-(n-2)} & z_{k-(n-3)} & \dots  & z_k
\end{bmatrix}$
$\mathbf{Q} = \mathbf{L} \mathbf{L}'$
$\mathbf{D} = (\sqrt{diag(\mathbf{Q})})^{-1}$
$\mathbf{R} = \mathbf{D} \mathbf{Q} \mathbf{D}'$
Here, I have guaranteed that $\mathbf{R}$ is a positive semi-definite correlation matrix that was generated from $\mathbf{z}$, which contains a bunch of randomly generated numbers. 
But I feel like there's a fundamental problem in the way I've set this up. I am using a numeric optimization procedure (in this specific case, python's scipy.optimize) to estimate $(n^2 - n)/2$ terms by searching within a $((n^2 - n)/2 +n)$-dimensional parameter space.
So there might be a bunch of places where the optimization surface is just flat because there are (almost?) infinite $\mathbf{z}$s that generate the same correlation matrix.
So finally, after all this, my question is: given the model stated up top, is there a better way to estimate all correlation terms through likelihood maximization in a way that guarantees positive semi-definiteness?
Any kind of guidance would be greatly appreciated! 
Edit
As suggested by @Mark L. Stone, Here is the implementation of the problem in Python with appropriate comments to make things a bit clearer, pointing out what parts of code are analogous to the formal/mathematical description I gave.
# Importing libraries used
import numpy as np
from scipy.stats import norm
from scipy.stats import multivariate_normal as mvn
from scipy.optimize import minimize

# Setting seed for replication
seed = 666
np.random.seed(seed)

# Number of dimensions in my jointly-distributed data. 
# Analogous to n.
ndim = 13

# Number of observations in the dataset.
# Analogous to m.
nobs = 1000

# Number of elements in parameter vector to be estimated.
# Analogous to k.
num_chol = int(((ndim*ndim)-ndim)/2+ndim)

# k-long vector of random numbers.
# Analogous to z.
true_chol_vec = norm(loc=0,scale=2.5).rvs(num_chol)

# Function that makes a covariance matrix using the random parameters.
def make_cov_mtx_from_chol_vec(chol_vec):
    chol_mtx = np.zeros((ndim,ndim))
    chol_mtx[np.tril_indices(ndim)] = chol_vec
    cov_mtx = np.dot(chol_mtx,chol_mtx.T)
    return(cov_mtx)

# Function that normalizes covariance matrix down to a correlation matrix.
def make_cor_mtx_from_cov_mtx(cov_mtx):
    stdevs = (1/np.sqrt(np.diag(cov_mtx))).reshape((ndim,1))
    cor_mtx = stdevs * cov_mtx * stdevs.T
    return(cor_mtx)

# Creating true covariance matrix. Analogous to Q. 
true_cov_mtx = make_cov_mtx_from_chol_vec(true_chol_vec)

# Creating true correlation matrix. Analogous to R.  
true_cor_mtx = make_cor_mtx_from_cov_mtx(true_cov_mtx)

# Mean of the jointly distributed data. Analogous to 0_n.
means = np.zeros(ndim)

# Generating correlated data to use in estimation. 
# Analogous to X.
joint_values = mvn.rvs(mean=np.zeros(ndim),cov=true_cor_mtx, size=nobs)

# Fixes cases where likelihoods are too small
prob_fix = 1e-10

# Log-likelihood function used in optimization. 
def neg_log_lik(params):
    cov_mtx_estim = make_cov_mtx_from_chol_vec(params)
    cor_mtx_estim = make_cor_mtx_from_cov_mtx(cov_mtx_estim)
    likelihood = mvn.pdf(x=joint_values,mean=means,cov=cor_mtx_estim, allow_singular=True)
    likelihood[likelihood < prob_fix] = prob_fix
    log_likelihood = np.log(likelihood)
    return(-log_likelihood.sum())

# Optimization starting values
start_params = norm(loc=0,scale=2.5).rvs(num_chol)

# Running optimization
param_optim = minimize(fun=neg_log_lik,
                       x0=start_params,
                       method="BFGS",
                       options={"maxiter":10000,
                                "disp":True})
print("\n\nTrue cor mtx:")
print(true_cor_mtx)

print("\n\nEstimated cor mtx:")
print(make_cor_mtx_from_cov_mtx(make_cov_mtx_from_chol_vec(param_optim["x"])))

In the implementation above, the estimation terminates after a single iteration and spits out a correlation matrix that isn't even close to what's expected. Is this because of the difference between the search-space dimension and the actual solution-space dimension?
In summary: what is the best way to estimate the correlation matrix of jointly correlated data using a numerical likelihood maximization approach? Is it possible that using this approach - where a $((n^2 - n)/2 + n)$-long parameter vector is used to fit $((n^2 - n)/2)$ correlation terms - might generate problems for the numerical optimization/search procedure? If so, how can that be avoided?
Thanks again!
 A: After almost giving up on this, I finally found a good paper that talks about unconstrained parametrizations for the covariance matrix.
The basic idea is to extract the upper triangle Cholesky decomposition and to apply several sequential spherical-coordinate transformations on each of the Cholesky's matrix's columns. It gets really messy and writing a generalized code is pretty tricky, but it's totally doable.
When performing the spherical transforms, you need to calculate two "sets" of parameters: $r$s and $\phi$s. Since you're working with a correlation matrix, however, you know that all of the $r$s are equal to 1, thus reducing the number of values for your optimization process.
Important Links:
Unconstrained parametrizations for variance-covariance matrices (link 2): gives you the general idea of using the Cholesky decomposition with spherical coordinate transformations. The step about spherical coordinates, however, isn't very clear here.
Parameterizing correlations: a geometric interpretation: Another look at a similar question.
n-dimension Spherical coordinates and the volumes of the n-ball in $\mathbb{R}^n$: How to generalize the spherical-coordinate transformation for n-dimensions.
Spherical Parameterization of Variance-Covariance Matrix in Mixed-Effects Regression: Question about the first paper.
Spherical Parametrization of a Cholesky Decomposition: Question about the first paper, specifically about the spherical coordinate transformation.
A: As R is a correlation matrix, you can fix the diagonal terms of the Cholesky. Q will not be a correlation matrix but R will be.
$m = (n^2-n)/2$
$\mathbf{z} =(z_1, z_2, ..., z_m)=m \in \mathbb{R}^m $
$\mathbf{L} = \begin{bmatrix}
    1    & 0      & 0      & \dots  & 0 \\
    z_1    & 1    & 0      & \dots  & 0 \\
    z_2    & z_3    & 1    & \dots  & 0 \\
    \vdots & \vdots & \vdots & \ddots & \vdots \\
    z_{m-(n-1)} & z_{m-(n-2)} & z_{m-(n-3)} & \dots  & 1
\end{bmatrix}$
$\mathbf{Q} = \mathbf{L} \mathbf{L}'$
$\mathbf{D} = (\sqrt{diag(\mathbf{Q})})^{-1}$
$\mathbf{R} = \mathbf{D} \mathbf{Q} \mathbf{D}'$
Some code in R to validate:
for (ind in 1:1000) {
  n <- sample (3:10,1)
  L <- diag(n)
  L[lower.tri(L, diag=FALSE)] <- rnorm((n^2 - n)/2)
  Q <- crossprod(L)
  D <- 1/sqrt(diag(Q))
  if (any(eigen(diag(D) %*% Q %*% diag(D))$values < 0)) stop("Error")
}

