在R中避免使用多个for循环计算矩阵

Question

在为回答一个地图问题生成一些虚假数据的过程中，我发现自己写出了以下内容：

# Generate some fake data
lat <- seq(-90, 90, by = 5)
lon <- seq(-180, 180, by = 10)
phi <- matrix(0, nrow = length(lat), ncol = length(lon))
i <- 1
for (l1 in lat) {
    j <- 1
    for (l2 in lon) {
        phi[i, j] <- (sin(pi * l1 / 180) * cos(pi * l2 / 180))^2
        j <- j+1
    }
    i <- i+1
}
phi <- 1500*phi + 4500  # scale it properly

很明显，那两个中心for循环不像我想要的R语言风格那样。看起来我应该能够使用 mapply 或其他工具来完成这项工作，但遗憾的是它将返回一个列表，而不是我想要的结果。其他 apply 函数似乎也不能做到我需要的事情。

我在这里错过了什么？

- Mike Wise

Answer 1

你应该尝试使用矩阵代数。无需使用任何来自apply函数族的函数：

lat <- seq(-90, 90, by = 5)
lon <- seq(-180, 180, by = 10)
1500 * tcrossprod(sin(pi * lat / 180), cos(pi * lon / 180))^2 + 4500

- Raad

---

原文在这里：http://stackoverflow.com/questions/35592266/raster-image-on-world-map-in-ggplot/35598289#35598289 - 我给你署名了。 - Mike Wise

干杯！必须得说下面的速度测试相当令人惊讶。 - Raad

除了对这个问题的兴趣外，这确实让我感到震惊。 - Mike Wise

虽然我仍然非常喜欢这个答案，并在我的问题中使用它，但我现在注意到下面的“outer”解决方案实际上更通用。 - Mike Wise

Answer 2

你可以使用outer。

   x = outer(lat, lon, FUN = function(x,y) {(sin(pi * x/180) * cos(pi * y /180))^2})
    identical(x * 1500 + 4500, phi)
# [1] TRUE

NBATrends的答案似乎比其他解决方案更快。以下是一些基准测试结果。

library(microbenchmark) 
microbenchmark(within(df, {
  phi <- (sin(pi * lat / 180) * cos(pi * lon / 180))^2
  phi <- 1500*phi + 4500
}), 1500 * tcrossprod(sin(pi * lat / 180), cos(pi * lon / 180))^2 + 4500, outer(lat, lon, FUN = function(x,y) {(sin(pi * x/180) * cos(pi * y /180))^2}),
((as.matrix(l1)%*%t(as.matrix(l2)))^2) * 1500 + 4500)
Unit: microseconds
                                                                                              expr     min       lq      mean   median       uq     max neval
 within(df, {     phi <- (sin(pi * lat/180) * cos(pi * lon/180))^2     phi <- 1500 * phi + 4500 }) 255.670 262.0095 270.50948 266.6880 277.7060 385.467   100
                                  1500 * tcrossprod(sin(pi * lat/180), cos(pi * lon/180))^2 + 4500  11.471  12.3770  22.30177  12.9805  13.5850 868.130   100
               outer(lat, lon, FUN = function(x, y) {     (sin(pi * x/180) * cos(pi * y/180))^2 }) 137.645 139.7590 144.39520 141.5700 145.1925 179.905   100
                                            ((as.matrix(l1) %*% t(as.matrix(l2)))^2) * 1500 + 4500  16.301  17.6595  20.20390  19.6215  20.5270  80.294   100

- Mamoun Benghezal

---

太多了。有趣。 - Mike Wise

现在我有时间仔细考虑，使用outer的答案在许多方面上是更好和更通用的答案，因为我们可以为x和y放置任意函数，并且crossprod实际上只执行乘积运算的函数。但是，在适用的情况下，crossprod速度要快得多，对于我的特定问题，crossprod非常适合，所以我不会调整正确的答案，并将其保留在这个注释中。 - Mike Wise

Answer 3

对于您的应用程序来说，线性代数可能更简单，因为您只需要逐元素相乘两个向量，这可以通过 v * u^T 完成。在 R 中，矩阵乘法是 %*%。

lat <- seq(-90, 90, by = 5)
lon <- seq(-180, 180, by = 10)

l1 <- sin(pi * lat / 180) 
l2 <- s(pi * lon/ 180)

# compute the matrix
phi <- as.matrix(l1)%*%t(as.matrix(l2))
# square each element of the matrix
phi <- phi^2
# scale properly
# square each element of the matrix
phi <- 1500*phi + 4500  

- Mathieu B

Answer 4

为什么要沉溺于矩阵结构和使用apply，当你可以向量化呢？

df <- expand.grid(lat = seq(-90, 90, by = 5),
                 lon = seq(-180, 180, by = 10))
df <- within(df, {
  phi <- (sin(pi * lat / 180) * cos(pi * lon / 180))^2
  phi <- 1500*phi + 4500
  })

您可以随时按照此处的说明转换回去。

- sebastian-c

Answer 5

我使用了sapply()，但我更喜欢outer()的解决方案：

#using sapply
phi_1 <- 
  t(
    sapply(lat, function(l1)
      sapply(lon, function(l2)(sin(pi * l1 / 180) * cos(pi * l2 / 180))^2))
  ) * 1500 + 4500

#compare result
identical(phi_1, phi)
# [1] TRUE

- zx8754