# Passing a cholesky decomposition for a matrix with constrained variances to an objective function

I am trying to optimize an objective function $$L(\theta)$$ in which some parameters that I aim to recover belong to a covariance matrix, $$\Sigma$$. $$\Sigma$$ has a unique structure, which includes ones on the diagonals: $$\Sigma = \begin{pmatrix} 1 & & \\ \sigma_{21} & 1 & \\ \sigma_{31} & \sigma_{32} & 1\\ \end{pmatrix} \\$$

If $$\Sigma$$ had no such restrictions, one would typically include the lower-triangular Cholesky decomposition, $$L$$, in $$\theta$$ and pass it to the objective function, and then recreate $$\Sigma = L L'$$ inside the objective function. This ensures that, regardless of the values that the optimizer tries, you always have a valid covariance matrix (symmetric and positive semidefinite).

In my case, because $$\Sigma$$ has ones on the diagonal, only three parameters are identified, which means that I can only pass three parameters from the Cholesky decomposition to my objective function. Let the Cholesky decomposition be written as $$L = \begin{pmatrix} x_{11} & 0 & 0\\ x_{21} & x_{22} & 0 \\ x_{31} & x_{32} & x_{33} \\ \end{pmatrix}$$

One option is to pass the parameters the off-diagonal parameters of $$L$$, $$x_{21}, x_{31}, x_{32}$$, and then set $$x_{11} = 1$$, $$x_{22} = \sqrt{1 - x_{21}^2}$$, and $$x_{33} = \sqrt{1 - x_{31}^2 - x_{32}^2}$$ inside the objective function. If I did this, then $$LL' = \Sigma.$$ However, this is suboptimal from an optimization standpoint, because $$x_{21}, x_{31}, x_{32}$$ cannot vary freely over the parameter space, but rather must satisfy a series of nonlinear constraints $$(x_{21}^2 < 1, x_{31}^2 + x_{32}^2 < 1)$$.

What is the best way of passing three parameters to the likelihood function, such that I can reconstruct $$\Sigma$$ inside the objective function without running into problems arising from taking the square root of a negative number?

• Are you using an environment with automatic differentiation (e.g. tensorflow, pytorch, jax) or are you planning on computing gradients manually? Commented Feb 1 at 22:23
• But yeah, generically, one option would be to parameterize $\mathbf{L}$ as a lower triangular matrix per usual, then map it onto the correlation matrices in two steps. 1) $\tilde{\Sigma}=\mathbf{L}\mathbf{L}^\top$, as usual, followed by 2) $\Sigma = \mathbf{D}\tilde{\Sigma}\mathbf{D}$, where $\mathbf{D}$ is a diagonal matrix with elements given by the inverse square root of the diagonal elements of $\tilde{\Sigma}$. Commented Feb 1 at 22:26
• I will not be using automatic differentiation. I'll either approximate the gradient numerically or compute it analytically. And thank you, your solution works. In the example above, I would pass $x_{21}, x_{31}, x_{32}$ to the objective function and set these equal to the lower triangular elements of L and set the diagonals equal to 1. Then, steps 1) and 2) above will provide a valid covariance matrix with the restrictions I want regardless of the values of $x_{21}, x_{31}, x_{32}$. Commented Feb 2 at 1:44
• cool I ask because it seems like it would be a pain to compute the gradient analytically, but if you do really only have a 3 by 3 matrix numerical gradient computation may suffice. Commented Feb 2 at 14:59
• since my comment seems to have answered the question I went ahead and made it an actual answer as CV ettiquette demands. Commented Feb 2 at 15:52

We can use parameters organized in the usual manner:

$$L = \begin{pmatrix} x_{11} & 0 & 0\\ x_{21} & x_{22} & 0 \\ x_{31} & x_{32} & x_{33} \\ \end{pmatrix}$$

And then map this onto the set of correlation matrices in a two step process.

First, mapping it onto the set of covariance matrices, per usual: $$\tilde{\Sigma} = \mathbf{L}\mathbf{L}^\top$$

Second, mapping this covariance matrix onto the set of correlation matrices by scaling the rows and columns by the inverse square root of the variances. Linear algebraically, this looks like forming a diagonal matrix $$\mathbf{D}$$: $$\textrm{diag}(\mathbf{D}) = \Bigg[ \frac{1}{\sqrt{\tilde\Sigma_{1,1}}}, \ldots, \frac{1}{\sqrt{\tilde\Sigma_{N,N}}}\Bigg]$$

and then pre and postmultiplying $$\tilde{\Sigma}$$ by this matrix:

$$\Sigma = \mathbf{D}\tilde\Sigma\mathbf{D}$$

When implementing this on a computer, it will be more efficient to simply scale the columns and rows of $$\tilde{\Sigma}$$, rather than explicitly computing the matrix multiplication.