Solving for most likely intersection of multiple multivariate Gaussians

Question

Is there a performant way to directly solve for the most likely intersection point (X, Y) of several multivariable Gaussians?

I've seen a few posts here that have asked how to solve for the intersection between two Gaussians - the concept is familiar to me. Right now it's not obvious to me aside from iterating and solving for two distributions at a time.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

mus = [np.array([[0.3],[0.7]]), 
       np.array([[0.3],[0.2]]), 
       np.array([[1.5],[0.6]])]        

covs = [np.array([[0.85, 0.3], [0.3, 0.25]]), 
        np.array([[0.7, -0.41], [-0.41, 0.25]]), 
        np.array([[0.5, 0.15], [0.15, 0.15]])]

cmaps = ["Reds", "Blues", "Greens"]

for m, cov, c in zip(mus, covs, cmaps):

    cov_inv = np.linalg.inv(cov) 
    cov_det = np.linalg.det(cov)

    x = np.linspace(-3, 3)
    y = np.linspace(-3, 3)
    X,Y = np.meshgrid(x,y)
    coe = 1.0 / ((2 * np.pi)**2 * cov_det)**0.5
    Z = coe * np.e ** (-0.5 * (cov_inv[0,0]*(X-m[0])**2 + (cov_inv[0,1] + cov_inv[1,0])*(X-m[0])*(Y-m[1]) + cov_inv[1,1]*(Y-m[1])**2))
    plt.contour(X,Y,Z, cmap = c)

Answer 1

You can do a LOT better than iterating between 2 solutions at a time. Realize that at every (x, y) point, you have a Z value for all 3 curves, and at the 3-way intersecting point, they are all equal (or within tolerance). And at other points, if you take the lowest Z of the curves, and move towards the center (mu_x, mu_y) of that curve, you are moving in an improving direction.

The below is an iterative algorithm that does that. There is certainly some meat on the bone in terms of possible enhancements. Notably, you could incorporate a "tolerance" for stopping conditions easily, or do some weighted average of the 2 lower z values instead of just the lowest to get the movement vector, or tinker with a larger step size.

Anyhow, this converges very rapidly for many different test starting points.

Code:

import numpy as np
import matplotlib.pyplot as plt

class Curve:
    # a convenience so we can avoid recomputations
    def __init__(self, mu, cov_inv, cov_det):
        self.mu = mu
        self.cov_inv = cov_inv
        self.cov_det = cov_det
        self.coe = 1.0 / ((2 * np.pi)**2 * cov_det)**0.5

    def z(self, x, y):
        Z = self.coe * np.e ** (-0.5 * (self.cov_inv[0,0]*(x-self.mu[0])**2 + \
            (self.cov_inv[0,1] + self.cov_inv[1,0])*(x-self.mu[0])*(y-self.mu[1]) + self.cov_inv[1,1]*(y-self.mu[1])**2))
        return Z

mus = [np.array([[0.3],[0.7]]), 
       np.array([[0.3],[0.2]]), 
       np.array([[1.5],[0.6]])]        

covs = [np.array([[0.85, 0.3], [0.3, 0.25]]), 
        np.array([[0.7, -0.41], [-0.41, 0.25]]), 
        np.array([[0.5, 0.15], [0.15, 0.15]])]

cmaps = ["Reds", "Blues", "Greens"]

curves = []
for m, cov, c in zip(mus, covs, cmaps):

    cov_inv = np.linalg.inv(cov) 
    cov_det = np.linalg.det(cov)

    x = np.linspace(-3, 3)
    y = np.linspace(-3, 3)
    X,Y = np.meshgrid(x,y)
    coe = 1.0 / ((2 * np.pi)**2 * cov_det)**0.5
    Z = coe * np.e ** (-0.5 * (cov_inv[0,0]*(X-m[0])**2 + (cov_inv[0,1] + cov_inv[1,0])*(X-m[0])*(Y-m[1]) + cov_inv[1,1]*(Y-m[1])**2))
    plt.contour(X,Y,Z, cmap = c)
    curves.append(Curve(m, cov_inv, cov_det))

# iterative algorithm...

pos = np.array((-1,2))
step_size = 0.1
num_steps = 100
footprints = [pos,]

for step in range(num_steps):
    zs = [ (curves[i].z(*pos), i) for i in range(len(curves))]
    zs.sort()  # sort by z value, lowest will be first
    c = curves[zs[0][1]]  # the curve to move toward
    vec = c.mu.T - pos
    move_vec = vec * (step_size/np.linalg.norm(vec))
    print(f'move: {move_vec} towards curve {zs[0][1]}')
    pos = pos + move_vec
    pos = pos.flatten()
    # check to see if we have backtracked, if so, shorten the step
    if len(footprints) > 1 and np.linalg.norm(pos - footprints[-2]) < step_size:
        #print(f'norm: {np.linalg.norm(pos-footprints[-2])}')
        step_size *= 0.5
    footprints.append(pos)

plt.plot([t[0] for t in footprints], [t[1] for t in footprints], c='k', lw=2)
plt.show()

Plot: