I codded my PDF function for the multivariate gaussian (3D) as such:

    def Gaussian3DPDF(v, mu, sigma):
      N = len(v)
      G = 1 / ( (2 * np.pi)**(N/2) * np.linalg.norm(sigma)**0.5 )
      G *= np.exp( - (0.5 * (v - mu).T.dot(np.linalg.inv(sigma)).dot( (v - mu) ) ) )
      return G

It works pretty well and it's quite fast. But I now want to fit my data to this function using a least square optimizer, and to be efficient I need to be able to pass a Mx3 matrix ( [[x0,y0,z0],[x1,...],...] )

However, my dot product will break down because now it's not 3 dot 3x3 dot 3 but 3xM dot 3x3 dot Mx3

I'm not very good with linear algebra.

Is there a trick I am not aware of ?

Thanks a lot

PS: Doing a for loop over each coordinate work but it is way way too slow for fitting on large number of data I have.

PPS: I found out about the scipy stats multivariate gaussian function. It works fine ! Though I'm still interested if anyone knows the answer ! :D


The code to do this in python without linear algebra:

#cube = np.array of dimention NxMxO

def Gaussian3DPDFFunc(X, mu0, mu1, mu2, s0, s1, s2, A):
    mu = np.array([mu0, mu1, mu2])
    Sigma = np.array([[s0, 0, 0], [0, s1, 0], [0, 0, s2]])
    res = multivariate_normal.pdf(X, mean=mu, cov=Sigma)
    res *= A
    res += 100
    return res

def FitGaussian(cube):
    # prepare the data for curvefit
    X = []
    Y = []
    for i in range(cube.shape[0]):
        for j in range(cube.shape[1]):
            for k in range(cube.shape[2]):
    bounds = [[3,3,3,3,3,3,50], [cube.shape[0] - 3, cube.shape[1] - 3, cube.shape[2] - 3, 30, 30, 30, 100000]]
    p0 = [cube.shape[0]/2, cube.shape[1]/2, cube.shape[2]/2, 10, 10, 10, 100]
    popt, pcov = curve_fit(Gaussian3DPDFFunc, X, Y, p0, bounds=bounds)
    mu = [popt[0], popt[1], popt[2]]
    sigma = [[popt[3], 0, 0], [0, popt[4], 0], [0, 0, popt[5]]]
    A = popt[6]
    res = multivariate_normal.pdf(X, mean=mu, cov=sigma)
    return mu, sigma, A, res

For linear algebra look at bellow ! Really Cool


1 Answer 1


Least squares optimizer has an elegant solution using linear algebra. You are solving the system $A\hat x=\hat b$, where be is A is your matrix ( [[1,x0,z0],[1,x1,y2],...] ), $b$ is a column of [z0; z1; ;..] and $x$ is a vector containing the estimated parameters which your solving for. The vector $b$ is NOT in the column space of $A$, so there is no solution, so you need to decompose the vector $b$ vector into the sum of the projection of $b$ onto the column space of the matrix $A$ and the orthogonal component $e$ given by the following:

\begin{equation} \label{2} b=proj_{Col(A)} + e \end{equation}

Where $e$ is a vector containing errors orthogonal to the column space of $A$.

Instead of solving $A x= b$, we solve the equation that best estimates $ b$.

\begin{equation} Ax=proj_{Col(A)} \end{equation}

Since, the $proj_{Col(A)}$ (read as the projection of b onto the column space of $A$) is in the column space of $A$ there will now be a solution to the system, where there wasn't one previously one! To find a the "best fit" solution start by combining the previous equations:

\begin{equation} A x= b - e \end{equation}

Here comes the trick! Multiply each term by $A^T$, which is the transposed matrix of A where the columns now become rows. \begin{equation} A^T A x=A^T b -A^T e \end{equation}

Where $e$ is orthogonal to the row space of $A^T$, and therefore $ e$ is in the null space of $A^T$. This means term $A^T e $ becomes the zero vector $ 0$. What's left is the least squares solution to $A x=b$ given by :

\begin{equation} A^T A x=A^T b \end{equation}

Now you want to know how to code this...

I will solve the 2 variable case:

Multiplying everything out we get:

Multiplying everything out we get: Multiplying everything out we get:

To solve for the estimators, the matrix should be augmented and row reduced.

The row reduction starts by switching row 1 and row 2. Then multiply row 1 by $-\frac{n}{\sum_{i=1}^{n} x_i}$ and add to row 2. This will result in a $0$ in the second row and first column. A total of two pivots for two rows means the matrix has full rank and $\hat b_0$ and $\hat b_1$ can be solved for.

  • $\begingroup$ Xqua, my question is off topic, but is this multivariate PDF Gaussian related to signal processing in flow Cytometers. The reason I ask is because your first two questions on this cite were about Gaussian fitting of a flat top voltage signal that I have seen last week in testing for saturation on flow Cytometers :) $\endgroup$
    – hwhorf
    Commented Jun 6, 2018 at 1:26
  • $\begingroup$ Whoa really cool answer !! Thanks ! It's so elegant !! @hwhorf Ahahah no this time, it's to find gaussian in a 3D microscopy volume ;) $\endgroup$
    – Xqua
    Commented Jun 6, 2018 at 2:09
  • $\begingroup$ Forming and solving the Normal Equations is not the best way to compute least squares stats.stackexchange.com/questions/343069/… . For some better ways, see stats.stackexchange.com/questions/160179/… . $\endgroup$ Commented Jun 6, 2018 at 2:36
  • $\begingroup$ yes, SVD is better for solving least squares than deriving equations. $\endgroup$
    – hwhorf
    Commented Jun 6, 2018 at 2:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.