是否有一个 R 函数来获取矩阵的相邻值的“加权”总和（对于不同的半径）？答案

【问题标题】：Is there an R function to get the "weighted" sum of neighboring values of a matrix (for varying radius)?是否有一个 R 函数来获取矩阵的相邻值的“加权”总和（对于不同的半径）？
【发布时间】：2019-03-24 05:23:13
【问题描述】：

我有一个矩阵 M（比如 r 行和 c 列），我想根据每个矩阵元素的邻居获得“加权”总和并创建一个新矩阵 M2。邻居一词可能在半径 1 内（在元胞自动机理论中通常称为摩尔邻域），或者半径可能与 1 不同，例如 2、3 等。

对于矩阵 M 中的特定单元格，请说在中间的某处。假设位置（i，j）；那么第（i，j）个单元格有“八个”邻居，

(i-1, j-1), (i-1, j), (i-1, j+1), (i, j-1), (i, j+1), (i+1, j-1), (i+1, j), (i+1, j+1).

我想创建一个矩阵 M2 来计算第 (i,j) 个单元格加上其八个相邻单元格的“加权”和。加权是基于单元之间的欧几里得距离进行的。比如，

exp(-sqrt(2))*M[i-1,j-1] + exp(-1)*M[i-1,j] + exp(-sqrt(2))*M[i-1,j+1] + exp(-1)*M[i,j-1] + M[i,j] + exp(-1)*M[i,j+1] + exp(-sqrt(2))*M[i+1,j-1] + exp(-1)*M[i+1,j] + exp(-sqrt(2))*M[i+1,j+1]

对所有单元重复相同的想法（沿边界的单元需要特殊处理，因为它们不一定有八个相邻单元）。上面的想法适用于半径 1，但我正在尝试开发的代码需要对任何半径都是通用的。

r <- 4
c <- 4

n <- r*c

(M <- matrix(1:n, r, c))

addresses <- expand.grid(x = 1:r, y = 1:c)

#Got this code in the same forum

z <- rbind(c(-1,0,1,-1,1,-1,0,1),c(-1,-1,-1,0,0,1,1,1))

get.neighbors <- function(rw) {
  # Convert to absolute addresses 
  z2 <- t(z + unlist(rw))
  # Choose those with indices within mat 
  b.good <- rowSums(z2 > 0)==2  &  z2[,1] <= nrow(M)  &  z2[,2] <= ncol(M)
  M[z2[b.good,]]
}

apply(addresses,1 , get.neighbors) # Returns a list with neighbors

M

本质上，半径 = 1 的 M2 必须是当前单元格加上相邻单元格的“加权”总和。当前当前单元格的权重始终为 1。

M = [ 1  5   9  13
      2  6  10  14
      3  7  11  15
      4  8  12  16]

M2 = [ 5.033  13.803 .... ....
       ....   ....   .... ....
       ....   ....   .... ....
       ....   ....   .... ....]

如何在 R 中获取矩阵 M2？如果半径大于1怎么办？我希望加权发生在两个 for 循环内，因此我可以在关闭两个 for 循环的代码中进一步使用计算得出的 [i,j] 单元格的加权和。

【问题讨论】：

请不要永远在问题中发布的代码中包含rm(list=ls())。这类似于建议某人运行您的代码，其中包含旧的format C:\ DOS 命令。
谢谢 r2evans。我会记住以后永远不要包含rm(list=ls())。

标签： r for-loop matrix nearest-neighbor adjacency-matrix

【解决方案1】：

我认为以下是您想要的加权和。我会以与@r2evans 类似的方式找到邻居。

wtd_nbrs_sum <- function(input_matrix,
                         radius,
                         weight_matrix)
{
  temp_1 <- matrix(data = 0,
                   nrow = nrow(x = input_matrix),
                   ncol = radius)
  temp_2 <- matrix(data = 0,
                   nrow = radius,
                   ncol = ((2 * radius) + ncol(x = input_matrix)))
  input_matrix_modified <- rbind(temp_2,
                                 cbind(temp_1, input_matrix, temp_1),
                                 temp_2)
  output_matrix <- matrix(nrow = nrow(x = input_matrix),
                          ncol = ncol(x = input_matrix))
  for(i in seq_len(length.out = nrow(x = input_matrix)))
  {
    for(j in seq_len(length.out = nrow(x = input_matrix)))
    {
      row_min <- (radius + (i - radius))
      row_max <- (radius + (i + radius))
      column_min <- (radius + (j - radius))
      column_max <- (radius + (j + radius))
      neighbours <- input_matrix_modified[(row_min:row_max), (column_min:column_max)]
      weighted_sum <- sum(neighbours * weight_matrix)
      output_matrix[i, j] <- weighted_sum
    }
  }
  return(output_matrix)
}

r <- 4
c <- 4
n <- r*c
M <- matrix(data = 1:n,
            nrow = r,
            ncol = c)
R <- 1
wts <- matrix(data = c(exp(x = -sqrt(x = 2)), exp(x = -1), exp(x = -sqrt(x = 2)), exp(x = -1), 1, exp(x = -1), exp(x = -sqrt(x = 2)), exp(x = -1), exp(x = -sqrt(x = 2))),
              nrow = 3,
              ncol = 3)

wtd_nbrs_sum(input_matrix = M,
             radius = R,
             weight_matrix = wts)
#>           [,1]     [,2]     [,3]     [,4]
#> [1,]  5.033856 13.80347 24.16296 23.89239
#> [2,]  8.596195 20.66391 34.43985 32.84175
#> [3,] 11.186067 24.10789 37.88383 35.43163
#> [4,]  9.748491 19.86486 30.22435 28.60703

^{由reprex package (v0.2.1) 于 2019 年 3 月 24 日创建}

希望这会有所帮助。

【讨论】：

谢谢@yarnabrina 我在更改R = 2 时收到错误Error in neighbours * weight_matrix : non-conformable arrays。我看到wts 矩阵对于R = 1 是3 x 3，但是随着R 的增加wts 的加权距离计算也会增加。非常感谢任何帮助。
@AshokKrishnamurthy 如果半径为 R，则需要将权重矩阵更新为 (2R+1)x(2R+1) 矩阵。所以，如果你尝试R=2，得到一个5x5 矩阵wts，然后运行wtd_nbrs_sum(input_matrix = M, radius = R, weight_matrix = wts)。

【解决方案2】：

编辑以包含加权和。

这样做可能有一个非常巧妙的技巧，但最直接（且可维护）的方法可能是一个简单的两个for-loop 实现。

M1 <- matrix(1:16, nr=4)
M1
#      [,1] [,2] [,3] [,4]
# [1,]    1    5    9   13
# [2,]    2    6   10   14
# [3,]    3    7   11   15
# [4,]    4    8   12   16

代码：

get_neighbors <- function(M, radius = 1) {
  M2 <- M
  M2[] <- 0
  nr <- nrow(M)
  nc <- ncol(M)
  eg <- expand.grid((-radius):radius, (-radius):radius)
  eg$wt <- exp(-sqrt(abs(eg[,1]) + abs(eg[,2])))
  for (R in seq_len(nr)) {
    for (C in seq_len(nc)) {
      ind <- cbind(R + eg[,1], C + eg[,2], eg[,3])
      ind <- ind[ 0 < ind[,1] & ind[,1] <= nr &
                    0 < ind[,2] & ind[,2] <= nc,, drop = FALSE ]
      M2[R,C] <- sum(M[ind[,1:2, drop=FALSE]] * ind[,3])
    }
  }
  M2
}

get_neighbors(M1, 1)
#           [,1]     [,2]     [,3]     [,4]
# [1,]  5.033856 13.80347 24.16296 23.89239
# [2,]  8.596195 20.66391 34.43985 32.84175
# [3,] 11.186067 24.10789 37.88383 35.43163
# [4,]  9.748491 19.86486 30.22435 28.60703

同样的东西，半径为 2：

get_neighbors(M1, 2)
#          [,1]     [,2]     [,3]     [,4]
# [1,] 12.44761 25.64963 31.73247 32.70974
# [2,] 18.57765 35.96237 43.33862 43.51911
# [3,] 20.09836 37.80643 45.18268 45.03982
# [4,] 17.51314 31.88500 37.96784 37.77527

还有一个简单的测试，如果使用半径 0，那么 M1 和 M2 应该是相同的（它们是）。

注意：这通常在基础 R 中执行很好，没有花哨地使用 apply 或其表亲。由于这是一个非常直接的启发式方法，因此可以很容易地使用 Rcpp 来实现它，以显着加快速度。

【讨论】：