Calculating a 3D gradient with unevenly spaced points

Question

I currently have a volume spanned by a few million every unevenly spaced particles and each particle has an attribute (potential, for those who are curious) that I want to calculate the local force (acceleration) for.

np.gradient only works with evenly spaced data and I looked here: Second order gradient in numpy where interpolation is necessary but I could not find a 3D spline implementation in Numpy.

Some code that will produce representative data:

import numpy as np    
from scipy.spatial import cKDTree

x = np.random.uniform(-10, 10, 10000)
y = np.random.uniform(-10, 10, 10000)
z = np.random.uniform(-10, 10, 10000)
phi = np.random.uniform(-10**9, 0, 10000)

kdtree = cKDTree(np.c_[x,y,z])
_, index = kdtree.query([0,0,0], 32) #find 32 nearest particles to the origin
#find the gradient at (0,0,0) by considering the 32 nearest particles?  

(My problem is very similar to Function to compute 3D gradient with unevenly spaced sample locations but there seemed to have been no solution there, so I thought I'd ask again.)

Any help would be appreciated.

Answer 1

这是一个 Julia 实现，可以满足您的要求

using NearestNeighbors

n = 3;
k = 32; # for stability use  k > n*(n+3)/2

# Take a point near the center of cube
point = 0.5 + rand(n)*1e-3;
data = rand(n, 10^4);
kdtree = KDTree(data);
idxs, dists = knn(kdtree, point, k, true);

# Coords of the k-Nearest Neighbors
X = data[:,idxs];

# Least-squares recipe for coefficients
 C = point * ones(1,k); # central node
dX = X - C;  # diffs from central node
 G = dX' * dX;
 F =  G .* G;
 v = diag(G);
 N = pinv(G) * G;
 N = eye(N) - N;
 a =  N * pinv(F*N) * v;  # ...these are the coeffs

# Use a temperature distribution of  T = 25.4 * r^2
# whose analytical gradient is   gradT = 25.4 * 2*x
X2 = X .* X;
C2 = C .* C;
T  = 25.4 * n * mean(X2, 1)';
Tc = 25.4 * n * mean(C2, 1)'; # central node
dT = T - Tc;       # diffs from central node

y = dX * (a .* dT);   # Estimated gradient
g = 2 * 25.4 * point; # Analytical

# print results
@printf "Estimated  Grad  = %s\n" string(y')
@printf "Analytical Grad  = %s\n" string(g')
@printf "Relative Error   = %.8f\n" vecnorm(g-y)/vecnorm(g)

该方法的相对误差约为1%。以下是几次运行的结果...

Estimated  Grad  = [25.51670916224472 25.421038632006926 25.6711949674633]
Analytical Grad  = [25.41499027802736 25.44913042322385  25.448202594123806]
Relative Error   = 0.00559934

Estimated  Grad  = [25.310574056859014 25.549736360607493 25.368056350800604]
Analytical Grad  = [25.43200914200516  25.43243178887198  25.45061497749628]
Relative Error   = 0.00426558

更新
我不太了解Python，但这里有一个似乎有效的翻译

import numpy as np
from scipy.spatial import KDTree

n = 3;
k = 32;

# fill the cube with random points
data = np.random.rand(10000,n)
kdtree = KDTree(data)

# pick a point (at the center of the cube)
point = 0.5 * np.ones((1,n))

# Coords of k-Nearest Neighbors
dists, idxs = kdtree.query(point, k)
idxs = idxs[0]
X = data[idxs,:]

# Calculate coefficients
C = (np.dot(point.T, np.ones((1,k)))).T # central node
dX= X - C                    # diffs from central node
G = np.dot(dX, dX.T)
F = np.multiply(G, G)
v = np.diag(G);
N = np.dot(np.linalg.pinv(G), G)
N = np.eye(k) - N;
a = np.dot(np.dot(N, np.linalg.pinv(np.dot(F,N))), v)  # these are the coeffs

#  Temperature distribution is  T = 25.4 * r^2
X2 = np.multiply(X, X)
C2 = np.multiply(C, C)
T  = 25.4 * n * np.mean(X2, 1).T
Tc = 25.4 * n * np.mean(C2, 1).T # central node
dT = T - Tc;       # diffs from central node

# Analytical gradient ==>  gradT = 2*25.4* x
g = 2 * 25.4 * point;
print( "g[]: %s" % (g) )

# Estimated gradient
y = np.dot(dX.T, np.multiply(a, dT))
print( "y[]: %s,   Relative Error = %.8f" % (y, np.linalg.norm(g-y)/np.linalg.norm(g)) )

更新 #2
我想我可以使用格式化的 ASCII 而不是 LaTeX 编写一些易于理解的内容...

给定一组 M 个 n 维向量（称为 b_k），找到一组
`coeffs（称它们为 a_k），它产生身份的最佳估计
`矩阵和零向量
`
` 中号
` (1) 分钟 ||E - I||，其中 E = 总和 a_k b_k b_k
` a_k k=1
`
` 中号
` (2) 最小值 ||z - 0||，其中 z = a_k b_k 之和
` a_k k=1
`
`
`请注意，不需要基向量 {b_k}
`是标准化的、正交的，甚至是线性独立的。
`
`首先，定义以下数量：
`
` B ==> 列为 b_k 的矩阵
` G = B'.B ==> B 乘以 B 的转置
` F = G @G ==> @ 表示哈达玛积
` v = diag(G) ==> 由 G 的 diag 元素组成的向量
`
`上述最小化相当于这个线性约束问题
`
` 求解 Fa = v
` st Ga = 0
`
`令 {X} 表示 X 的 Moore-Penrose 逆。
`那么线性问题的解可以写成：
`
` N = I - {G}.G ==> 投影到 G 的零空间
` a = N 。{FN} 。v
`
`这些系数的用途是它们允许你写
`张量场导数的非常简单的表达式。
`
`
`设 D 为 del（或 nabla）运算符
`并且 d 是关于中心（又名第 0）节点的差分运算符，
`这样，对于任何标量/向量/张量 Y，我们有：
` dY = Y - Y_0
`
`设 x_k 为第 k 个节点的位置向量。
对于我们的基本向量，取
` b_k = dx_k = x_k - x_0。
`
`假设每个节点都有一个与其关联的字段值
`（例如温度），并假设二次模型 [关于 x = x_0]
` 对于字段 [g=梯度，H=hessian，“:”是双点积]
`
` Y = Y_0 + (x-x_0).g + (x-x_0)(x-x_0):H/2
` dY = dx.g + dxdx:H/2
` D2Y = I:H ==> Y 的拉普拉斯算子
`
`
`在第 k 个节点评估模型
`
` dY_k = dx_k.g + dx_k dx_k:H/2
`
`乘以 a_k 并求和
`
` 嗯
` 总和 a_k dY_k = 总和 a_k dx_k.g + 总和 a_k dx_k dx_k:H/2
` k=1 k=1 k=1
`
` = 0.g + I:H/2
` = D2Y / 2
`
“因此，我们有拉普拉斯算子的二阶估计
`
` 中号
` 圈(Y_0) = 总和 2 a_k dY_k
` k=1
`
`
`现在用线性模型玩同样的游戏
` dY_k = dx_k.g
`
`但是这次乘以 (a_k dx_k) 并求和
`
` MM
` 总和 a_k dx_k dY_k = 总和 a_k dx_k dx_k.g
` k=1 k=1
`
` = 免疫球蛋白
` = g
`
`
一般来说，中心节点的导数可以估计为
`
` 中号
` D#Y = 总和 a_k dx_k#dY_k
` k=1
`
` 中号
` D2Y = 总和 2 a_k dY_k
` k=1
`
` 哪里
` # 代表{点、叉或张量}积
` 产生 Y 的 {div、curl 或 grad}
` 和
` D2Y 代表 Y 的拉普拉斯算子
` D2Y = D.DY = 圈(Y)