如何在R语言中计算n维数据数组的梯度

我认为肯定有现成的R包函数，但始终找不到能对**数据点数组（非函数）**通过中心差分计算梯度的工具，要求能处理所有非边缘点，且支持任意n维数据集。

示例：3维网格点的梯度计算 #

下面是生成5x5x5网格数据点，并计算(3,3,3)点梯度的示例：

response_data <- array(NA, dim = c(5,5,5))
for(i in 1:dim(response_data)[1]) {
  for(j in 1:dim(response_data)[2]) {
    for(k in 1:dim(response_data)[3]) {
      response_data[i,j,k] <- i*j*k
    }
  }
}
  
x_values <- 1:5
y_values <- 1:5
z_values <- 1:5

# 计算(3,3,3)点的中心差分梯度
partial_x <- (response_data[4, 3, 3] - response_data[2, 3, 3])/
  (x_values[4] - x_values[2])
partial_y <- (response_data[3, 4, 3] - response_data[3, 2, 3])/
  (y_values[4] - y_values[2])
partial_z <- (response_data[3, 3, 4] - response_data[3, 3, 2])/
  (z_values[4] - z_values[2])

已有尝试与局限 #

• 编辑1：pracma::gradient可处理1维或2维数组，但不支持n维数组。
• 编辑2：我自己编写了3维数组的梯度计算代码，但想优化为支持任意n维的版本，目前卡在如何动态选择维度进行循环：

response_data <- array(NA, dim = c(5,5,5))
for(i in 1:dim(response_data)[1]) {
  for(j in 1:dim(response_data)[2]) {
    for(k in 1:dim(response_data)[3]) {
      response_data[i,j,k] <- i*j*k
    }
  }
}

# 将所有坐标存入列表
coord_vals <- list(x_values = 1:5, y_values = 1:5, z_values = 1:5)

# 创建空数组存储偏导数，每个元素对应一个自变量的偏导
pds <- rep(list(array(NA, dim = sapply(coord_vals, length))),
           length(coord_vals))
names(pds) <- names(coord_vals)

# 计算梯度：第一维度
dimension <- 1
for(j in 1:length(coord_vals[[2]])) {
  for(k in 1:length(coord_vals[[3]])) {
    pds[[dimension]][ , j, k] <- 
      pracma::gradient(F = response_data[, j, k],
                       h1 = coord_vals[[dimension]])
  }
}
# 第二维度
dimension <- 2
for(i in 1:length(coord_vals[[1]])) {
  for(k in 1:length(coord_vals[[3]])) {
    pds[[dimension]][i, , k] <- 
      pracma::gradient(F = response_data[i, , k],
                       h1 = coord_vals[[dimension]])
  }
}
# 第三维度
dimension <- 3
for(i in 1:length(coord_vals[[1]])) {
  for(j in 1:length(coord_vals[[2]])) {
    pds[[dimension]][i, j, ] <- 
      pracma::gradient(F = response_data[i, j, ],
                       h1 = coord_vals[[dimension]])
  }
}

内容的提问来源于stack exchange，提问作者mikeblazanin