使用共轭梯度法用 Numpy 数组在 3D 中拟合高斯有困难

Question

本质上，我将有一些 3d 强度分布，在我的问题的简化版本中，我将以 3D space 中的一个点为中心的一些高斯分布，由 3D numpy 数组表示，定义为：

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

我们的想法是将一些高斯分布拟合到一些未知的强度分布，改变它们的振幅，以有限数量的网格点为中心（选择强度高于最小阈值，在本例中为 0.5 的网格点），由于问题的凸性以及缩放问题，使用梯度下降算法。

我们要最小化，其中x_v表示所有的网格点，而G表示在给定网格点[=]处以网格点x_i为中心的高斯值34=] 振幅 s_i。这具有以下梯度：

而这个渐变是用下面的代码实现的：

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

然而，对于作为相应网格点处强度测量值的振幅初始估计，此分析梯度不同于使用有限差分方案发现的梯度，因此 cg 算法不起作用：

但是调试了几天，一直找不到问题的根源。为不同问题计算梯度的类似方法用于许多相同的代码，并且此实现在 2D 中也适用，没有 z 轴贡献。一定有一些基本的东西我错过了，但我不确定是什么。

我在 https://pastebin.com/rhs4tasZ 上有一份完整的测试代码副本，以供需要的任何其他信息，或供任何想运行自己使用此代码的人使用

Answer 1

问题已解决，我现在将关闭此问题。我在 synth_image

中提取了错误的 z 参数

使用共轭梯度法用 Numpy 数组在 3D 中拟合高斯有困难

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

python

arrays

numpy

gaussian

gradient-descent