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

# 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)

# 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

# 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)

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

# Running optimization
param_optim = minimize(fun=neg_log_lik,
print("\n\nTrue cor mtx:")

print("\n\nEstimated cor mtx:")

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!

  • $\begingroup$ You don't seem to have told us what optimization problem you are solving with scipy.optimize. Despite your long presentation, I have essentially no idea what you've done, other than relying on a Choelsky factorization to ensure positive semidefiniteness. You apparently populate a Choelsky factor with a random entries, but is that just an initialization (starting value) for numerical optimization via some unstated optimization problem formulation? $\endgroup$ Aug 5, 2018 at 17:32
  • $\begingroup$ Good point. I'll edit the original post to reflect that. Thx for the heads up! $\endgroup$
    – Felipe D.
    Aug 5, 2018 at 18:40
  • $\begingroup$ What output is displayed when the optimizer terminates after one iteration? By virtue of your approach (even if no mistakes), you may have introduced spurious saddle points. .Is there a reason why you can;t just form the empirical covariance matrix - do you have unclean, inconsistent data (not all variable components measured together)? If not, then other than roundoff errors, empirical covariance should be psd. You can adjust eigenvalues or use other methods to adjust an almost psd "covariance" or correlation matrix to be psd, or have minimum eigenvalue. $\endgroup$ Aug 5, 2018 at 20:00
  • $\begingroup$ The output I get is this: "Optimization terminated successfully. Current function value: 23025.850930 Iterations: 0 Function evaluations: 1277 Gradient evaluations: 1" $\endgroup$
    – Felipe D.
    Aug 5, 2018 at 20:16
  • $\begingroup$ The thing is that this is a small part in other larger estimation problems I have to solve. The main problem I deal with is the Generalized Ordered Probit Model, where I have multiple (discrete) ordered outcomes and the error terms are jointly distributed. In this problem, we estimate the influence of a bunch of exogenous covariates as well as the correlation terms between the errors (more info here). So I tried to translate the simplest version of the problem to a clean-cut context to present it here. $\endgroup$
    – Felipe D.
    Aug 5, 2018 at 20:18

3 Answers 3


There is I believe a much easier way to parameterize correlation matrices vis the cholesky decomposition. Form the matrix

\begin{bmatrix}1 & 0 & 0 & \ldots & 0 & 0 \\ x_{21} & 1 & 0 & \ldots & 0 & 0 \\ x_{31} & x_{32} & 1 & \ldots & 0 & 0 \\ \ldots \\ x_{n-1,1} & x_{n-1,2} & x_{n-1,3} & \ldots & 1 & 0 \\ x_{n1} & x_{n2} & x_{n3} & \ldots & x_{n,n-1} & 1 \\ \end{bmatrix}

Now normalize this by dividing each row considered as a vector by its norm. This results in the Cholesky decomposition of the correlation matrix. There are $n(n-1)/2$ $x$'s. Note that there are no constraints necessary on the $x$'s. To get the Cholesky decomposition of the covariance matrix multiply the $i$'th row by $\exp(y_i)$. So there are $n$ parameters for the $y$'s. Again note that there are no constraints necessary for the $y$'s. Notice that there is a natural starting value for the $x$'s, that is $x_i=0$ for all $i$. However the problem is singular at that point so you should start with small random values for the $x$'s.

I have used this parameterization for a lot of models and it seems to perform well. Consider a $4\times 4$ matrix so that the dimension of x is 4x3/2=6. Let

$$ x=(0.6379 , 0.4829 , 0.2525 ,-0.2435,-0.09326 ,-0.3671)$$

Then the matrix before normalizing is \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0.6379 & 1 & 0 & 0 \\ 0.4829 & 0.2525 & 1 & 0 \\ -0.2435 & -0.09326 & -0.3671 & 1 \\ \end{bmatrix}

and after normalizing it is

\begin{bmatrix} 1 & 0 & 0 & 0 \\ 0.5378 & 0.8431 & 0 & 0 \\ 0.4241 & 0.2217 & 0.8781 & 0 \\ -0.2221 & -0.08504 & -0.3347 & 0.9118 \\ \end{bmatrix}

If you multiply this matrix times its transpose you get the correlation matrix.

\begin{bmatrix} 1 & 0.5378 & 0.4241 & -0.2221 \\ 0.5378 & 1 & 0.4149 & -0.1911 \\ 0.4241 & 0.4149 & 1 & -0.4069 \\ -0.2221 & -0.1911 & -0.4069 & 1 \end{bmatrix}

  • $\begingroup$ Thanks for the answer, Dave! I'm trying this out with a small example and I'm running into some issues, which makes me think I didn't quite get your explanation. For example, consider a case with 3 dimenstions: $$ A_{lowmat}=\begin{bmatrix} 1 & 0 & 0 \\ 0.5 & 1 & 0 \\ 0.25 & 0.5 & 1 \end{bmatrix} $$ Here, when I divide each row by its norm, I get this: $$ \begin{bmatrix} 1 & 0 & 0 \\ 0.447 & 0.894 & 0 \\ 0.218 & 0.463 & 0.872 \end{bmatrix} $$ But multiplying this matrix by its transpose does not recover the original $A_{lowmat}$ matrix. Am I off? $\endgroup$
    – Felipe D.
    Apr 22, 2019 at 20:49
  • $\begingroup$ Also, you mention $y$ parameters, but it's not clear to me where they came from. Thanks! $\endgroup$
    – Felipe D.
    Apr 22, 2019 at 20:57
  • $\begingroup$ OK, if you are just doing the correlation, not the covariance matrix then there are no y's. $\endgroup$ Apr 22, 2019 at 23:29
  • $\begingroup$ Ah, ok. But this is a one-way process, right? You can only go from random draws (parametrized version) to an actual correlation matrix (unparametrized version), right? In other words, this method doesn't allow you to start off with a given correlation matrix and find the decomposition that, multiplied by its own transpose, would recover the original correlation matrix. Is that right? $\endgroup$
    – Felipe D.
    Apr 23, 2019 at 6:37
  • 1
    $\begingroup$ Well I'm not sure why you think that, but the answer is no. The map from the x's to the correlation matrices is a well behaved diffeomorphism, so 1 to 1. $\endgroup$ Apr 24, 2019 at 22:48

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.


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")

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.