加速倍数为14,000倍，也就是。获胜的计算机科学

在我作为科学软件开发人员的工作中，我倾向于编写大量代码。大多数没有读过计算机科学学位的人倾向于认为CS“只是”把代码扔到屏幕上，然后运行它。我与我的许多同事和具有其他背景的同事保持着良好的工作关系。…。物理、气候科学、生物学等等。但是，当谈到开发软件时，我得到的明显印象是，人们会想，“嘿，这能有多难呢？我们只需写下几条指令，告诉计算机我们想要做什么，然后点击执行按钮，然后‘砰！’一声，我们就得到了答案！“。

这种思路的问题在于，编写与您认为的意思不符的指令非常容易。例如，计算机可能完全无法解释您的程序。此外，在不实际执行程序的情况下，几乎无法判断程序是否会真正终止。而且有很多、很多、很多种方法来编写程序，使其执行起来“很慢”。“慢”是…。真的很慢。好像要花一辈子甚至更长的时间才能真正执行。最后一个问题是我在阅读没有接受过CS教育的人编写的软件时最常遇到的问题。解决这个问题是我的工作。

关于CS，人们没有意识到的一件事是，它教会了你关于计算的理论，可计算性(也就是说，我们真的能计算什么吗？我们经常理所当然地认为我们可以！)，算法复杂性，以及所有的知识、逻辑和分析技术，并帮助您编写一个将在最少的时间或使用最少的空间内运行的程序。

请允许我向您展示一个巨大的优化示例，这是我对一个同事编写的简单脚本所做的优化。

在气候科学方面，我们做了大量的缩减工作。我们从粗尺度的全球气候模式网格中获取温度和降水读数，并将其映射到细尺度的局部网格。假设全局网格为50x25，本地网格为1000x500。对于局部网格中的每个网格单元，我们希望知道它对应于全局网格中的哪个网格单元。

考虑这个问题的一个简单方法是，我们想要最小化L[n]和G[n]之间的距离。因此，执行搜索的简单方法是：

对于每个局部像元L[i]：对于每个全局像元G[j]：计算L[i]和G[j]之间的距离查找集合L[i]*G中的最小距离返回最小值的索引。

这看起来很简单。然而，如果你仔细观察，你会注意到你必须做很多额外的工作。根据输入的大小来查看算法。

对于每个局部像元L[i]：#对每个全局像元G[j]执行此操作L次：#执行此操作L x G次计算L[i]和G[j]之间的距离(D)#执行此L x G次查找集合中的最小距离d[i*j]#读取G个像元L次(开销L x G)查找其像元与最小值匹配的索引#读取G个像元L次(开销L x G)。

Obs.lon<；-ncvar_get(nc.obs，'；lon&39；)obs.lat<；-ncvar_get(nc.obs，'；lat<；)n.lon<；-length(obs.lon)n.lat<；-length(obs.lat)obs.lat<；-Matrix(obs.lat，nrow=n.lon，nol=n.lat，byrow=true)obs.lon<；-Matrix(obs.lon，nrow=n.lon，ncol=n.lat)obs.time<；-netcdf.endar(nc.obs)gcm.lon<；-ncvar_get(nc.gcm，'；lon&39；)-360gcm.lat<；-ncvar_get(nc.gcm，&39；lat&39；)gcm.lats<；-Matrix(gcm.lat，ncol=length(gcm.lat)，nrow=length(gcm.lon)，byrow=true)gcm.lons<；-Matrix(gcm.lon，ncol=length(gcm.lat))gcm.lons.lats<；-Matrix(gcm.lon，ncolumn=length(gcm.lat)，nrow=length(gcm.lat)，nrow=length(gcm.lon)，nrow=length(gcm.lon)，nrow=length(gcm.lat)，nrow=length(gcm.lat)。-cbind(c(gcm.lons)，c(gcm.lats))#找出哪些GCM网格框与每个精细网格点相关联#将搜索限制在10度。X 10°。Neighborhooddxy<；-10mdist<；-function(x，y)Apply(abs(Sweep(data.Matrix(Y)，2，data.Matrix(X)，'；-'；)，1，sum)nn<；-list()for(i in seq_long(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<；-this(gcm.lims)n.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数据，有1068x510个单元格(总共544,680个)，假设我们的GCM有50x25个单元格(总共1250个单元格)。因此，“某个计算单元”中内循环的成本为：

其中项是对应于计算两点之间的距离、找到最小点和找到数组索引的成本的常量。实际上，我们不太关心常量项，因为它们不受输入大小的影响。所以我们可以把它们聚集在一起，称之为成本；

这看起来很多，但真的吗？电脑很快，对吧？如果我们运行这个朴素的实现，情况大概是这样的：

最终耗时1668秒，略低于半小时。

>；source('；BCCA/naive.implementation.R'；)500 1000 1500 2000 2500 3000...543000 543500 544000 544500[1]&34；用户系统运行时间为1668.868 8.926 1681.728。

但是我们需要花30分钟吗？事情是这样的。我们正在比较两个网格，这两个网格都有大量我们没有利用过的结构。例如，粗格网和细格网中的纬度和经度都是按顺序排序的。因此，如果你想搜索一个数字，你不需要查看每一个数字。您可以使用平分算法，在该算法中，您可以查看中间的点，然后决定要搜索阵列的哪一半。那么搜索整个空间只需要花费搜索空间的日志(以2为基数)。

另一个我们没有利用的主要结构是纬度在维度中重复，而经度在维度中重复。所以不用做一次运算，我们可以做几次。这是一个巨大的优化。

对于每个LOCAL[x]：BICET_SEARCH(LOCAL[x]，GLOBAL[x])对于每个LOCAL[y]：BIVECT_SEARCH(LOCAL[y]，GLOBAL[y])返回每个维度的搜索结果的二维格网。

##对*排序后的向量v#返回最接近xfind最近的元素的数组索引。-function(x，v){if(length(V)==1){return(1)}if(length(V)==2){return(which.min(abs(v-x)}mid<；-ceuling(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(粗点，细点){return(sApply(fine.point，find.neest，粗点))}##取纬度和经度的精细比例尺(例如ANUSPLINE)栅格##并找到对应于粗略比例尺(例如GCM)栅格的索引##因为搜索实质上是2维中的最小化距离##我们实际上可以单独在每个维度中独立地搜索(其中##是巨大的优化，使运行时x+y而不是x*y)和##然后重构索引以创建完整的gridregrid.coarse.to.fin.lt；-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)##索引Xi<；二维网格-Matrix(xi，ncol=length(fine.lats)，nrow=length(fine.lons)，byrow=F)yi<；-Matrix(yi，ncolumn=length(fine.lons)，byrow=T)return(list(xi=xi，yi=yi)}。

每个平分搜索的成本是输入大小的对数。这次我们的输入大小被划分为X和Y空间，所以我们将为Global、Local、X和Y使用and。

加上我们的数字，我们估计成本为553,076英镑。55.3万听起来比6.8亿好多了。我们在运行时看到了这一点吗？

>；ptm<；-proc.time()；rv<；-regrid.coarse.to.Fine(gcm.lat，gcm.lon，obs.lat，obs.lon)；打印('；已用时间&39；)；打印(pro.time()-ptm)[1]&34；用户系统已用时间0.117 0.000 0.117>；STR(房车)列表2$Xi：Num[1：1068，1：510]15 15 15...$yi：Num[1：1068，1：510]13 13 13。

0.117秒。半小时前我们花了近半个小时，现在我们花了一秒钟多一点的时间。

Soooooo…。我知道我受过做这类工作的训练，我的工作就是知道如何做这类事情。但即使是我也对这种加速是如此显著感到惊讶和自以为是。这是14000倍的加速。

此脚本过去需要很长时间，因此必须将其输出保存到磁盘，并由科学家手动检查，然后才能继续。现在你一眨眼的功夫就可以计算出来。这是一个我们必须做几百次的计算，这为我们节省了几天到几周的计算时间。它增加了与系统交互的能力，帮助我们从科学家的Time…中获得更多价值。他们不会无所事事地等待计算完成。它就这么做了。

我应该强调的是，这些史诗般的性能改进是在没有购买任何更大的计算机系统、没有并行化或增加复杂性的情况下出现的…。事实上，更快算法的代码实际上更简单、更可重用！只要阅读代码并对计算复杂性有一定的了解，这几乎是一个全面的胜利。

由Disqus提供支持的博客评论