Difficulty with fitting Gaussians in 3D with Numpy arrays using conjugate gradient method

Question

Essentially, I will have some 3d intensity distribution, where in a simplified version of my problem I will have some Gaussian centered at a point in 3D space, represented by a 3D numpy array, defined as:

I = np.zeros((ny, nx, nz))
tolerance = 1e-4 # minimum value of Gaussian for compact representation
sigma = x[1]-x[0] # Gaussian width
max_dist = sigma*(-2*np.log(tolerance))
di = np.ceil(max_dist/dx) # maximum distance in compact representation, in index format

# Create intensity field/true Gaussian
# this exists separately as its own function synth_I() where [0] is instead for each particle [i]
ix = round((xp[0] - x[0]) / dx) # index of particle center
iy = round((yp[0] - y[0]) / dy)
iz = round((yp[0] - y[0]) / dz)
iix = np.arange(max(0, ix - di), min(nx, ix + di), 1, dtype=int) # grid points with nonzero intensity values
iiy = np.arange(max(0, iy - di), min(ny, iy + di), 1, dtype=int)
iiz = np.arange(max(0, iz - di), min(nz, iz + di), 1, dtype=int)
ddx = dx * iix - xp[0] # distance between particle center and grid point
ddy = dy * iiy - yp[0]
ddz = dz * iiz - zp[0]
gx = np.exp(-1 / (2 * sigma ** 2) * ddx ** 2) # 1D Gaussian
gy = np.exp(-1 / (2 * sigma ** 2) * ddy ** 2)
gz = np.exp(-1 / (2 * sigma ** 2) * ddz ** 2)
gx = gx[np.newaxis,:, np.newaxis]
gy = gy[:,np.newaxis, np.newaxis]
gz = gz[np.newaxis, np.newaxis, :]
I[np.ix_(iiy, iix, iiz)] = I[np.ix_(iiy, iix, iiz)] + gy*gx*gz

The idea is to fit some a series of Gaussians, varying their amplitudes, centered at a finite number of grid points (chosen as those with an intensity above a minimum threshold, 0.5 in this case) to some unknown intensity distribution, using a gradient descent algorithm due to the convexity of the problem but also due to scaling concerns.

We want to minimize where x_v represents all the grid points, while G represents the value of a Gaussian centered at grid point x_i at a given grid point x_v with amplitude s_i. This has the following gradient:

And this gradient is implemented with the following code:

def diff_and_grad(s):  # ping test this:
    part_params = np.concatenate((xpart, ypart, zpart, s)) # for own code
    # create an intensity field as a combination of Gaussians as above
    synthI = synth_I_field_compact(part_params, nd, sigma, x, y, z)
    Idiff = I - synthI # difference in measurements
    f = 0.5 * np.sum(np.power(Idiff, 2)) # objective function
    g = np.zeros(Np) # gradient
    for i in range(0, Np):
        ix = round((xpart[i] - x[0]) / dx)
        iy = round((ypart[i] - y[0]) / dy)
        iz = round((zpart[i] - z[0]) / dz)
        iix = np.arange(max(0, ix - di), min(nx, ix + di), 1, dtype=int)
        iiy = np.arange(max(0, iy - di), min(ny, iy + di), 1, dtype=int)
        iiz = np.arange(max(0, iz - di), min(nz, iz + di), 1, dtype=int)
        ddx = dx * iix - xpart[i]
        ddy = dy * iiy - ypart[i]
        ddz = dz * iiz - zpart[i]
        gx = np.exp(-1 / (2 * sigma ** 2) * ddx ** 2)
        gy = np.exp(-1 / (2 * sigma ** 2) * ddy ** 2)
        gz = np.exp(-1 / (2 * sigma ** 2) * ddz ** 2)
        Id = Idiff[np.ix_(iiy, iix, iiz)]
        g[i] = -Id.dot(gz).dot(gx).dot(gy) # gradient is -product of local intensity difference with gaussian centered at grid point
    return f, g

However, for an initial estimate of amplitude taken as the value of the intensity measurements at the corresponding grid points, this analytical gradient differs from one found with a finite difference scheme, and thus the cg algorithm does not work:

However, after a couple days of debugging, I have been unable to find the source of the problem. A similar method of calculating the gradients for a different problem is used with much of the same code, and as well this implementation works in 2D without the z axis contributions as well. There must be something fundamental I am missing, but I am not sure what.

I have a copy of the full test code at https://pastebin.com/rhs4tasZ for any additional information needed, or for anyone who wants to run this code themselves

Answer 1

Fixed the problem, I'll close this question now. I extracted the wrong parameters for z in synth_image