For each $j = 1,\dots,N$, let $\mu_j \in \mathbb{R}^N$ denote a known column vector, $\Sigma_j \in \mathbb{R}^{N\times N}$ a known covariance matrix, and $\theta_j \in \mathbb{R}$ an unknown parameter, which we may write as a vector $\theta = (\theta_1, \dots, \theta_N)$. Let $\mathcal{N}(0,1)$ be a multivariate normal distribution of dimension $N$ and let $X \in \mathbb{R}^Nm$ be a multivariate normal distribution defined by
\begin{align} X &= (\sqrt{\theta_1} \mathcal{N}(0, \Sigma_1) + \theta_1 \mu_1) + \dots + (\sqrt{\theta_N} \mathcal{N}(0, \Sigma_N) + \theta_N \mu_N)\\ &= \mathcal{N}(\theta_1 \mu_1, \theta_1 \Sigma_1) + \dots + \mathcal{N}(\theta_N \mu_N, \theta_N \Sigma_N) \\ &= \mathcal{N}\left( \sum_{j=1}^N \theta_j\mu_j, \sum_{j=1}^N \theta_j \Sigma_j \right) \,. \end{align}
Given only one sample $X_s = x = (x_1, \dots, x_N) \in \mathbb{R}^N$ drawn from $X$, I am trying to find a closed form solution for the maximum likelihood estimate $\hat{\theta}$, whilst being aware that this may also simply be intractable, not to mention non-identifiable.
Given the various forms $X$ as written above, I have attempted finding $\frac{\partial f}{\partial \theta_j}$ of the likelihood function
\begin{align} f(x|\theta) = \frac{1}{(2\pi)^{N/2} \mathrm{det}(\sum_{j=1}^N \theta_j \Sigma_j)} \exp\left( -\frac{1}{2}(x- \sum_{j=1}^N \theta_j\mu_j)^T (\sum_{j=1}^N \theta_j \Sigma_j)^{-1} (x- \sum_{j=1}^N \theta_j\mu_j) \right) \end{align} and whilst this derivative was possible to calculate I couldn't see a way to manipulate it to solve for $\theta_j$ in any sense.
I also investigated using the first form of $X$ and using convolution identities to perhaps simplify the derivative, but this ultimately lead to solving something of the form $(g(x,\theta) f_1(x,\theta)) * \dots * f_N(x,\theta) = 0$, which was not much fun to deal with either.
Letting $P = [\mu_1, \dots, \mu_N] \in \mathbb{R}^{N \times N}$, a simple method of moments gives an estimate $\tilde{\theta} = P^{-1}x$ (where $x$ is a column vector), though this is to be taken with a grain of salt as this doesn't encapsulate the fact that $\theta$ is also a parameter of the covariance. It is not clear to me if there is a way to reconcile this. If the closed form solution fails to be fruitful then I may resort to testing how numerical solutions to the MLE compare to that moments estimate.
Any solutions, tips, insights or references are greatly appreciated!