I have a discrete representation of a surface/height-map $z = f(x,y)$ that i want to model as a mixture of Gaussians (please take probability distributions out of your mind for a moment). More precisely,

$ z = f(x) = \sum\limits_{i=1}^{n_g} a_i \; g(x; c_i, \Sigma_i)$


  • $x \in R^d$
  • $g(x; c_i, \Sigma_i) = exp\{-\frac{(x-c_i)^T\Sigma_i^{-1}(x-c_i)}{2} \}$,
  • $a_i \geq 0$ is the max-height of the $i$-th Gaussian component,
  • $c_i \in R^d$ is the center of the $i$-th Gaussian component, and
  • $\Sigma_i$ is the covariance of the $i$-th Gaussian component that determines its support in space and its orientation.

Given the data of the surface, I want to fit this model, i.e estimate the parameters $[...,a_i, c_i, \Sigma_i, ...]$ given some reasonable initial guesses by minimizing the least squares objective function given below:

$E(a,c,\Sigma) = \sum\limits_{j=1}^N \left(z_j - f(x_j; a, c, \Sigma)\right)^2$


  • $a_i \geq 0$,
  • $\Sigma_i$ should be symmetric ($\Sigma_i = \Sigma_i^T$), and positive-definite ($u^T\Sigma_i u > 0 \; \forall{u} \in R^d$).

I would like to emphasize that I want to optimize for the full covariance matrices $\Sigma_i$ (even the diagonal won't suffice).

Approach so far

I will describe below the approach I have been venturing on so far but I would really appreciate any other ideas in fitting this model to data or at least point me to the right literature to read (Googling takes me to so much literature on fitting a mixture model to a probability distribution using EM). I might have a huge number of components in the mixture. So it will be nice if you could supplement your answers with comments on computational complexity, the quality of the solution (how far from global minimum) and sensitivity to initialization (how close should the initial guesses be).

My initial idea is to use iterative minimization algorithms such as the Gauss-Newton algorithm or the Levenberg–Marquardt algorithm.

So I started off by deriving the derivatives of the error/objective function $E(a,c,\Sigma)$ with respect to each of its parameters:

$\begin{eqnarray} \frac{\partial E}{\partial a_i} & = & -2\sum\limits_{j=1}^N \left(z_j - f(x_j; a, c, \Sigma)\right) \frac{\partial f(x_j; a, c, \Sigma)}{\partial a_i} \\ & = & -2\sum\limits_{j=1}^N \left(z_j - f(x_j; a, c, \Sigma)\right) g(x_j; c_i, \Sigma_i) \end{eqnarray}$

$\begin{eqnarray} \frac{\partial E}{\partial c_i} & = & -2\sum\limits_{j=1}^N \left(z_j - f(x_j; a, c, \Sigma)\right) \frac{\partial f(x_j; a, c, \Sigma)}{\partial c_i} \\ & = & -2\sum\limits_{j=1}^N \left(z_j - f(x_j; a, c, \Sigma)\right) \; a_i \frac{\partial g(x_j; c_i, \Sigma_i)}{\partial c_i} \\ & = & -2 \sum\limits_{j=1}^N \left(z_j - f(x_j; a, c, \Sigma)\right) \; a_i \; g(x_j; c_i, \Sigma_i) \; \Sigma_i^{-1} (x_j - c_i) \end{eqnarray}$

$\begin{eqnarray} \frac{\partial E}{\partial \Sigma_i} & = & -2\sum\limits_{j=1}^N \left(z_j - f(x_j; a, c, \Sigma)\right) \frac{\partial f(x_j; a, c, \Sigma)}{\partial \Sigma_i} \\ & = & -2\sum\limits_{j=1}^N \left(z_j - f(x_j; a, c, \Sigma)\right) \; a_i \frac{\partial g(x_j; c_i, \Sigma_i)}{\partial \Sigma_i} \\ & = & -2 \sum\limits_{j=1}^N \left(z_j - f(x_j; a, c, \Sigma)\right) \; a_i \; g(x_j; c_i, \Sigma_i) \; \frac{\Sigma_i^{-T} (x_j - c_i)(x_j - c_i)^T\Sigma_i^{-T}}{2} \end{eqnarray}$

Now while I minimize the error function $E(a,c,\Sigma)$ using one of the iterative minimization algorithms, I want to enforce that the covariance matrices $\Sigma_i$ are symmetric and positive definite. And I don't know how to enforce that constraint? (should I try to optimize for elements on one side of the diagonal which is sort of projecting it to the domain of symmetric matrices but not necessarily positive definite?) I am guessing this is an already solved problem, at least reasonably. It would be great if you could point me to the literature for this problem so I can move forward.

Alternatively, coming from a computer-vision + machine-learning background, I can't resist but to think about applying the Expectation-Maximization algorithm. Intuitively, if I knew which Gaussian component each $(x,y,z)$ location belongs to (in a soft sense) then I can compute the parameters of the Gaussian components. And, if I know the parameters of the Gaussian mixture, I can compute the degree to which each $(x,y,z)$ location belongs each component in the mixture. And, Expectation-Maximization algorithms in the end boils down to minimization of an objective function using coordinate descent. So it would be nice if you could point me to any literature you can point me to that takes the EM approach to this problem.

  • $\begingroup$ Why do you say "please take probability distributions out of your mind for a moment" ? Why fit the height $a_i$? If you model a mixture of gaussians, the max height of the $i$th component of the mixture distribution will be the normalization constant. Your $a_i= p_i/ (i\text{th normalization constant})$ where $p_i$ is the proportion of observations belonging to the $i$th component. See the gaussian mixture example here that can be easily extended to your problem. $\endgroup$
    – jsk
    May 23 '14 at 7:48
  • $\begingroup$ @jsk my problem is different. My data is not a bunch of (x,y,z) points where in you are trying to apply GMM clustering using the EM algorithm. Let me simplify to 1D so you will understand. Think of a curve generated by summing up of gaussians each with a different center, different sigmas, and different height at the center. Given such a curve and given the number of gaussians, i want to solve the inverse problem of estimating the parameters of the each of the gaussian components in the mixture. Now, imagine it in higher dimensions (>=1). $\endgroup$
    – cdeepakroy
    May 23 '14 at 15:13
  • $\begingroup$ You may ask what changes in higher dimensions. In higher dimensions, each of the gaussian components will have a covariance matrix which tells you how it is oriented (given by the principal axes or eigen vectors of the covariance matrix) and its support in space (given by the eigen values of the covariance matrix which correspond to variances along the principal axes). $\endgroup$
    – cdeepakroy
    May 23 '14 at 17:24
  • $\begingroup$ I don't see the difficulty in setting the constraint: In the end, the matrix multiplication can be understood as a sum of products and you can just reuse the model parameters of Sigma - you simply have fewer parameters. However, how to find the right number of Gaussians? $\endgroup$
    – Make42
    May 18 '18 at 10:40

Your idea is exactly what I would recommend. To enforce the positive definiteness of the covariance matrices, just use the Cholesky decomposition $\Sigma_i = L_iL_i^T$, and minimise the objective with respect to $L_i$ instead of $\Sigma_i$. I refer to this article for more details. Also, you may try to optimise the log-likelihood instead of the mean square error.

  • $\begingroup$ that looks like an interesting idea and i feel good i wasnt off track. But for cholesky decomposition $L_i$ should be a lower traingular matrix with real and positive diagonal entries. How do i enforce that constraint? What comes to my mind is to put in bounds constraints on the elements of the matrix $L_i$. May be set the lower and upper bounds of the upper diagonal elements both to zero so the optimizer doesnt tamper with them. And, set the lower bound of the diagonal elements to zero and upper bound to infinity. Does this make sense? $\endgroup$
    – cdeepakroy
    May 23 '14 at 23:26
  • $\begingroup$ Also, which optimizer can handle bound constraints? If you know, can you point me to a reference which optimizes an objective function with symmetric and positive definite constraint enforced. That would make my mind at rest? $\endgroup$
    – cdeepakroy
    May 23 '14 at 23:29
  • $\begingroup$ I have updated my answer to provide a reference. $\endgroup$
    – Emile
    May 24 '14 at 9:58

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.