Estimating correlation matrix using numeric likelihood maximization

Question

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!

Answer 1

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}

Answer 2

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.

Answer 3

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