14,000倍加速或计算机科学胜利

作为科学软件开发人员，我做了很多编程工作。而且其他科学领域的大多数人倾向于认为编程只是“抛出”代码并运行代码。我与许多同事有着良好的工作关系，包括来自其他国家的同事……物理学，气候学，生物学等。但是，在软件开发方面，您会得到他们独特的印象，他们认为：“嘿，可能是困难吗？我们只需写下一些有关计算机应该执行的操作的说明，然后按“运行”按钮，我们就完成了-我们得到了答案！”



问题是，编写不代表您的想法的指令非常容易。例如，程序可能完全违反计算机的解释。除了，从字面上看，没有办法断定一个程序是否不执行就终止。而且有很多很多方法可以大大降低程序的执行速度。我的意思是...真的慢了。放慢速度，这样您将需要一生或更多的时间才能完成。这通常是由未经计算机教育的人（即来自其他领域的科学家）编写的程序发生的。我的工作是修复此类程序。



人们不理解计算机科学会教给您计算理论，算法的复杂性和可计算性（也就是说，我们实际上可以计算出什么吗？我们常常认为它可以做到！）将在最短的时间内或最少的资源使用中执行的代码。



让我向您展示由我的一位同事编写的一个简单脚本中的巨大优化示例。



我们在气候学上有很多研究。我们从大规模的全球气候模型中获取温度和降水读数，并将其与精细的本地网格进行比较。假设全局网格为50x25，本地网格为1000x500。对于本地网格中的每个网格单元，我们想知道它对应于全局网格中的哪个网格单元。



表示问题的一种简单方法是最小化L [n]和G [n]之间的距离。事实证明这样的搜索：

    L[i]:
      G[j]:
        L[i]  G[j]
       L[i] * G
    

似乎很简单。但是，如果您仔细观察，我们会做很多额外的工作。从输入的大小来看算法。

    L[i]:                           #  L 
      G[j]:                        #  L x G 
        L[i]  G[j]                 #  L x G 
       d[i*j]              #  G  L  (L x G)
   ,       #  G  L  (L x G)

代码看起来像这样：

obs.lon <- ncvar_get(nc.obs, 'lon')
obs.lat <- ncvar_get(nc.obs, 'lat')
n.lon <- length(obs.lon)
n.lat <- length(obs.lat)

obs.lats <- matrix(obs.lat, nrow=n.lon, ncol=n.lat, byrow=TRUE)
obs.lons <- matrix(obs.lon, nrow=n.lon, ncol=n.lat)
obs.time <- netcdf.calendar(nc.obs)

gcm.lon <- ncvar_get(nc.gcm, 'lon')-360
gcm.lat <- ncvar_get(nc.gcm, 'lat')
gcm.lats <- matrix(gcm.lat, ncol=length(gcm.lat), nrow=length(gcm.lon),
                   byrow=TRUE)
gcm.lons <- matrix(gcm.lon, ncol=length(gcm.lat), nrow=length(gcm.lon))
gcm.lons.lats <- cbind(c(gcm.lons), c(gcm.lats))

# Figure out which GCM grid boxes are associated with each fine-scale grid point
# Confine search to 10 deg. x 10 deg. neighbourhood

dxy <- 10
mdist <- function(x, y)
    apply(abs(sweep(data.matrix(y), 2, data.matrix(x), '-')), 1, sum)
nn <- list()
for (i in seq_along(obs.lons)) {
    if((i %% 500)==0) cat(i, '')
    gcm.lims <- ((gcm.lons.lats[,1] >= (obs.lons[i]-dxy)) &
                 (gcm.lons.lats[,1] <= (obs.lons[i]+dxy))) &
                ((gcm.lons.lats[,2] >= (obs.lats[i]-dxy)) &
                 (gcm.lons.lats[,2] <= (obs.lats[i]+dxy)))
    gcm.lims <- which(gcm.lims)
    nn.min <- which.min(mdist(c(obs.lons[i], obs.lats[i]),
                        gcm.lons.lats[gcm.lims,]))
    nn[[i]] <- gcm.lims[nn.min]
}
nn <- unlist(nn)

它看起来像一个简单的算法。计算距离然后找到最小值很容易。但是在这种实现方式中，随着局部像元数量的增加，我们的计算成本随着其乘积与全局网格中像元数量的乘积而增加。加拿大的ANUSPLIN数据包含1,068×510个单元格（总计544,680个）。假设我们的GCM包含50x25的单元（总共1250个）。因此，以“某些计算单位”为单位的内部循环的成本为：

$$

$$(c_0 \cdot L \times G) + (c_1 \cdot L \times G) + (c_2 \cdot L \times G)$$

成员在哪里 $$$c$是对应于计算两点之间的距离，找到最小点并找到数组索引的成本的常数。实际上，我们并不在乎常量成员，因为它们不依赖于输入的大小。因此，您可以将它们加起来并估算计算成本；

$$

$$(c \cdot L \times G)$$

因此，对于这组输入，我们的成本为 $$$544,680 \times 1,250 = 680,850,000$...



6.8亿。



似乎很多……虽然，谁知道呢？电脑速度很快，对不对？如果我们运行一个幼稚的实现，实际上它将比1668秒快一点，这比半小时少一点。

> source('BCCA/naive.implementation.R')
500 1000 1500 2000 2500 3000 ... 543000 543500 544000 544500 [1] "Elapsed Time"
    user   system  elapsed 
1668.868    8.926 1681.728 

但是程序应该运行30分钟吗？那就是问题所在。我们正在比较两个网格，每个网格都有大量我们从未使用过的结构。例如，两个网格中的纬度和经度均按排序顺序。因此，找到一个数字，您不需要遍历整个数组。您可以使用减半算法-查看中间的点，并确定我们要数组的哪一半。然后搜索整个空间将以整个搜索空间的底两个对数为准。



我们没有使用的另一个重要结构是纬度在维度上$$$x$和经度-在维度中 $$$y$... 因此，代替执行操作$$$x \times y$ 因为我们可以做到 $$$x + y$时间。这是一个巨大的优化。



伪代码看起来像什么？

  local[x]:
    bisect_search(local[x], Global[x])

  local[y]:
    bisect_search(local[y], Global[y])

 2d      

在代码中：

## Perform a binary search on the *sorted* vector v
## Return the array index of the element closest to x
find.nearest <- function(x, v) {
    if (length(v) == 1) {
        return(1)
    }
    if (length(v) == 2) {
        return(which.min(abs(v - x)))
    }
    mid <- ceiling(length(v) / 2)
    if (x == v[mid]) {
        return(mid)
    } else if (x < v[mid]) {
        return(find.nearest(x, v[1:mid]))
    }
    else {
        return((mid - 1) + find.nearest(x, v[mid:length(v)]))
    }
}

regrid.one.dim <- function(coarse.points, fine.points) {
    return(sapply(fine.points, find.nearest, coarse.points))
}

## Take a fine scale (e.g. ANUSPLINE) grid of latitudes and longitudes
## and find the indicies that correspond to a coarse scale (e.g. a GCM) grid
## Since the search is essentially a minimizing distance in 2 dimensions
## We can actually search independently in each dimensions separately (which
## is a huge optimization, making the run time x + y instead of x * y) and
## then reconstruct the indices to create a full grid
regrid.coarse.to.fine <- function(coarse.lats, coarse.lons, fine.lats, fine.lons) {
    xi <- regrid.one.dim(gcm.lon, obs.lon)
    yi <- regrid.one.dim(gcm.lat, obs.lat)
    ## Two dimensional grid of indices
    xi <- matrix(xi, ncol=length(fine.lats), nrow=length(fine.lons), byrow=F)
    yi <- matrix(yi, ncol=length(fine.lats), nrow=length(fine.lons), byrow=T)
    return(list(xi=xi, yi=yi))
}

每次平分搜索的成本是输入大小的对数。这次，我们将输入大小分为X和Y空间，因此我们将使用$$$G_x, G_y, L_x$ 和 $$$L_y$ 用于全局，本地，X和Y。

$$

$$cost = L_x \times log_2 G_x + L_y \times log_2 G_y + L_x \times L_y$$

成本为553076，好吧，553k的声音比680m好得多。这将如何影响运行时间？

> ptm <- proc.time(); rv <- regrid.coarse.to.fine(gcm.lat, gcm.lon, obs.lat, obs.lon); print('Elapsed Time'); print(proc.time() - ptm)[1] "Elapsed Time"
   user  system elapsed 
  0.117   0.000   0.117 
> str(rv)
List of 2
 $ xi: num [1:1068, 1:510] 15 15 15 15 15 15 15 15 15 15 ...
 $ yi: num [1:1068, 1:510] 13 13 13 13 13 13 13 13 13 13 ...
> 

0.117秒。过去需要将近半小时的时间现在仅需十分之一秒。

> 1668.868 / .117
[1] 14263.83

Soooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo，但即使我来说，惊喜和令人印象深刻的加速度是多少。约1.4万次。



以前，脚本太慢了，以至于将结果保存到磁盘上，以供科学家进行手动检查，然后再继续。现在，一切都在眨眼间就计算完了。我们进行了数百次这样的计算，因此最终节省了几天甚至几周的计算时间。而且我们有机会更好地与系统进行交互，从而使科学家的工作时间更具收益……他们不会闲着，等待计算结束。一切都准备就绪。



我必须强调，这些史诗般的性能改进不需要购买任何大型计算机系统，并行化或增加复杂性……实际上，更快的算法代码甚至更简单，更通用！只需通过仔细阅读代码并具有一定的计算复杂性知识，即可完全获得胜利。