Closed Form Solution for Gaussian Matrix which is Convex Combination?

Question

I already asked a pretty similar question here, but the answer was inconclusive and now this related problem has come up again here.

My problem is as follows, I have a $2n$-dimensional multivariate normal distribution with the following variance-covariance matrix: $$ \begin{bmatrix} \sigma_{x}^{2}[(\lambda)A_{n}+(1-\lambda)\mathbb{I}_{n}] & 0_{n x n} \\ 0_{n x n} & \sigma_{y}^{2}[(\lambda)A_{n}+(1-\lambda)\mathbb{I}_{n}] \\ \end{bmatrix}, $$ where $\lambda$ is unknown in $[0,1]$ [or if necessary, it can be made in $(0,1)$]; $\sigma_{x}^{2}$ and $\sigma_{y}^{2}$ are unknown variance factors greater than 0; $A_{n}$ is an $(n x n)$ known, positive semi-definite matrix with a diagonal of $1$s; $\mathbb{I}_{n}$ is the n-dimensional identity matrix.

Additionally, its mean is the concatenation of two n-length vectors of identical values $\mu_{x}$ and $\mu_{y}$.

Is there a closed form solution for any of the maximum likelihood estimates of $\lambda$, $\sigma_{x}^{2}$, $\sigma_{y}^{2}$, $\mu_{x}$, or $\mu_{y}$?

EDIT: Xi'an posted a comment saying this is unrelated to Gaussian mixtures, but isn't my problem fitting a linear combination of two Gaussian distributions? Here is a reformulation which I think highlights this: $$ (\lambda)\begin{bmatrix} \sigma_{x}^{2}A_{n} & 0_{n x n} \\ 0_{n x n} & \sigma_{y}^{2}A_{n} \\ \end{bmatrix}+ (1-\lambda)\begin{bmatrix} \sigma_{x}^{2}\mathbb{I}_{n} & 0_{n x n} \\ 0_{n x n} & \sigma_{y}^{2}\mathbb{I}_{n} \\ \end{bmatrix} $$

EDIT 2: mlofton asks whether my covariance is fixed or a function of t. I think it's probably either depending on the exact question. $A_{n}$ represents the expected correlation between observations at different points in space. Specifically, the correlation between $x_{i}$ and $x_{j}$ can be found at the entry in the ith row and jth column. These correlation values are function the spatial distance. But, we know the spatial distances and function which converts them into correlations, so these values are all fixed. I describe the way these distances are measured and converted into correlations here and here. But basically the variables evolve along a network of paths, and the covariance between two observations is proportional to the amount of shared path they have. So I'd say that because we know the correlations, their is no longer any explicit dependency on time, distance, or anything else. As for the $0_{nxn}$ matrices, in some versions of this problem they can be replaced by $\sigma^{2}c_{i}K_{n}$ where $c_{i}$ is a correlation factor between $X$ and $Y$, $K_{n}$ is a fully known matrix which describes the structure of the correlation, and $\sigma^{2}$ can be calculated from $\sigma_{x}^2$ and $\sigma_{y}^2$.

EDIT 3: I'd like to double check the formula I found online. For the unbiased variance $\sigma_{x}^2$ I should evaluate $$ \frac{1}{(n-1)} (x_{observed}-\mu_{x}1_{n})^{T}(\lambda A_{n}+(1-\lambda)\mathbb{I}_{n})^{-1}(x_{observed}-\mu_{x}1_{n}) $$ but probably I need the mle estimate instead right? So I simply replace $(n-1)$ with $(n)$ in the bottom part of the first fraction?

And then to calculate $\mu_{x}$, I suppose the formula is this? $$ \bigg(1_{n}^{T}\big((\lambda) A_{n}+(1-\lambda)\mathbb{I}_{n}\big)^{-1}1_{n}\bigg)^{-1}1_{n}^{T}\big((\lambda) A_{n}+(1-\lambda)\mathbb{I}_{n}\big)^{-1}x_{observed} $$

EDIT 4: I never verified the above formulas (in EDIT 3) because I found a black box solution for the mean and variance stuff. The dense grid search idea (proposed by Madden) is working really well though!

Answer 1

The answer to your question

Given $\lambda$, we can calculate $\mu_x, \mu_y, \sigma_x, \sigma_y$ using standard techniques, which you can find by searching "Weighted Least Squares" and adapting some of that math.

Unfortunately, it's not clear to me how optimization wrt $\lambda$ might be achieved in closed form. With a bit of work, we find that the log likelihood pertaining to $\lambda$ may be written as:

$$ - \Bigg[\sum_{i=1}^N \log\big(\lambda(\gamma_i-1)+1\big)\Bigg] - \frac{1}{2} \Bigg[ \sum_{i=1}^N \frac{\frac{u_{x,i}^2}{\sigma_x^2}+\frac{u_{y,i}^2}{\sigma_y^2}}{\lambda(\gamma_i-1)+1} \Bigg] \, , $$

where $\gamma_i$ is the $i$th eigenvalue of $\mathbf{A}_N$, $u_{x,i} = \mathbf{v}_i^\top(\mathbf{x}-\mu_x)$, and $u_{y,i} = \mathbf{v}_i^\top(\mathbf{y}-\mu_x)$, where $\mathbf{v}_i$ denotes the eigenvector of $\mathbf{A}_N$ corresponding to the $i$th eigenvalue.

To my knowledge, such an equation must be solved numerically. (Indeed, this is essentially the problem of determining a nugget parameter in Gaussian process regression, which is performed numerically in practice). So even if we knew everything else, we wouldn't be able to calculate $\lambda$ in closed form.

However, we can develop efficient numerics for this likelihood in practice. Given that you can compute every other parameter in closed form, it is simply a matter of varying $\lambda$ over a dense grid of $[0,1]$, computing the other parameters and the likelihood value at which the maximum is achieved, and taking that $\lambda$ which has the highest likelihood as our estimate. To do this efficiently, you will want to precompute the quantities I list above in order to use only one eigendecomposition of $\mathbf{A}_N$ to compute as many $\lambda$ evaluations as desired.

Why isn't this a normal mixture?

If we have two objects $a$ and $b$, and $\lambda\in[0,1]$, we call $\lambda a + (1-\lambda) b$ a convex combination.

What you have written is normal distribution with a covariance matrix which is a convex combination.

Gaussian mixtures also involve convex combinations, but it's the density itself which is a convex combination, not its covariance parameter.

So if we wanted to do a Gaussian mixture of the two normal densities $N(\mathbf{x}; \mu_x, \sigma^2_x \mathbf{A}_N)$ and $N(\mathbf{x}; \mu_x, \sigma^2_x \mathbf{I})$, that would yield the following density (up to proportionality):

$$ f(\mathbf{x}) \propto \lambda \big| \sigma^2_x \mathbf{A}_N \big| ^{-1} e^{-\frac{1}{2} (\mathbf{x}-\mu_x)^\top (\sigma^2_x \mathbf{A}_N)^{-1}(\mathbf{x}-\mu_x)} + (1-\lambda)\big| \sigma^2_x \mathbf{I} \big| ^{-1} e^{-\frac{1}{2} (\mathbf{x}-\mu_x)^\top (\sigma^2_x \mathbf{I})^{-1}(\mathbf{x}-\mu_x)} \, . $$

This is distinct from the density you have proposed, which (ignoring $\mathbf{y}$) is given by:

$$ f(\mathbf{x}) = \big| \sigma^2_x \big(\lambda\mathbf{A}_N+(1-\lambda)\mathbf{I}\big) \big| ^{-1} e^{-\frac{1}{2} (\mathbf{x}-\mu_x)^\top \big(\lambda\mathbf{A}_N+(1-\lambda)\mathbf{I}\big)^{-1}(\mathbf{x}-\mu_x)} \, . $$