栅格数据和矢量数据矢量化_从矢量化的诅咒中解放出数据科学家的思想-分页Julia进行救援...

栅格数据和矢量数据矢量化

如今，大多数数据科学家都使用Python或R作为主要编程语言。 在我今年早些时候遇到Julia之前，这也是我的情况。 Julia保证了与静态类型的编译语言（如C）相当的性能，同时保留了解释语言（如Python，R或Matlab）的快速开发功能。 此性能是通过即时（JIT）编译实现的。

Julia无需解释代码，而是在运行时编译代码。 尽管JIT编译已经存在了一段时间（例如Matlab在2002年引入了它 ），但是Julia的设计目的是考虑到JIT编译的性能。 类型稳定性和多次调度是Julia中的关键设计概念，这使其在竞争中脱颖而出。 如果您想了解更多信息，加利福尼亚大学的数据科学计划会提供一个非常不错的笔记本来解释这些概念。

在某个时候，我们开始使用解释语言来处理大型数据集（我猜数据集越来越大，我们一直使用相同的工具）。 我们了解到，为了提高性能，我们希望避免循环和递归。

相反，我们希望使用矢量化操作或专门的实现，这些操作或专业实现将数据结构（例如数组，数据帧）作为输入并在单个调用中处理它们。 我们这样做是因为在解释型语言中，每次执行一条指令时，我们都要支付开销。 当我很高兴用R进行编码时，它涉及到一套避免循环和递归的策略，很多时候，我们的努力都指向“如何避免解释型语言的陷阱？”。

我到了要编写C函数来解决R脚本瓶颈的地步，尽管性能得到明显改善，但使用R的优势却逐渐消失了。 那是我开始寻找替代品时发现的茱莉亚。

在本文中，我们将从解决R中的一个简单问题开始，我将尝试说明使用解释语言进行编程时的心态和局限性。 然后，我们将使用Julia来解决相同的问题，展示思维方式完全不同，以及如何才能实现类似C的性能。

矢量化途径

有时，不清楚如何使用矢量化获得最佳性能。

让我们考虑在给定点向量的情况下寻找点对的所有组合之间的距离的问题。 为了简单起见，点具有一维，我们将使用L1距离。 给定输入[5、3、9、1]，预期输出为[2、4、4、6、2、8]（由| 5–3 |，| 5–9 |，| 5–1 |产生， | 3-9 |，| 3-1 |和| 9-1 |）。

在R中使用stats包的dist函数可以解决此问题：

> as.vector(dist(c(5 , 3 , 9 , 1 ), method= "manhattan" ))
[ 1 ] 2 4 4 6 2 8

R的实现返回一个距离矩阵，我们将其转换为向量。

让我们假装dist不可用。 您将如何在R中进行编码？

用基于循环的方法解决这个问题很简单：我们需要一个外部循环来迭代该对中的第一个元素，并需要一个内部循环来第二个元素。

由于L1距离是对称的（| ab | = | ba | ），我们只需要对一半的组合进行此操作，同时避免计算点到它们自身的距离（矩阵的对角线）。 R中基于循环的实现如下所示：

distances_among_pairs_loop <-function (v) {
  nin = length(v)
  nout = (nin * nin - nin) %/% 2 #length of the output vector
  dists = vector(mode=typeof(v), length=nout) #preallocating output
  k = 1 #current position in the output vector
  for (i in seq_len(nin- 1 )) {
    a = v[i]
    for (b in v[seq(i+ 1 , nin)]) {
      dists[k] = abs(a-b)
      k = k + 1
    }
  }
  dists
}

如果您使用R进行编程，您会说这不是正确的R代码…循环很慢并且会产生非常冗长的代码……应该有一种向量化的方法。

丢弃循环的一种方法是通过生成具有所有对组合的向量。 基本上，我们正在找到距离矩阵下三角的（行，列）坐标。

然后，我们可以使用单个矢量化运算来计算距离。 我们知道我们要付出内存代价，但是我们希望向量化能带来回报。

distances_among_pairs_rep <-function (v) {
  row <- rep( 1 :length(v), each=length(v))  #e.g 1 1 1 2 2 2 3 3 3
  col <- rep( 1 :length(v), times=length(v)) #e.g 1 2 3 1 2 3 1 2 3
  lower_tri <- which(row > col) #e.g. (2,1) (3,1) (3,2)
  abs(v[row[lower_tri]] - v[col[lower_tri]])
}

另一种选择是使用R的outer函数生成一个矩阵，该矩阵具有两点的所有组合（包括冗余组合）之间的差。

然后，我们只需要检索矩阵下三角部分的绝对值。 我们可能不愿意，因为我们要进行的操作要多出大约2倍（并将它们存储在内存中），但这确实可以使代码更清晰，可读性更好。

distances_among_pairs_outer <-function (v) {
  dists <- outer(v, v, '-' )
  abs(dists[lower.tri(dists)])
}

当我们沿着向量化路径行进时，代码变得越来越紧凑，但是变得更快了吗？ 我们正在通过更多的内存和更多的操作来交换紧凑性……预测这些实现中最有效的实现并非易事。

当我们已经有一个简单的基于循环的解决方案（在编译语言中很难克服）时，我们可以继续尝试找出避免R陷阱的最佳方法是什么。

因此，如果此功能对性能至关重要，则可以在C，CPP或Fortran中实现它。 通过Rcpp库与R集成的原始实现的CPP转换如下所示：

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector distances_among_pairs_cpp(NumericVector v) {
    int nin = v.size();
    int nout = (nin * nin - nin) / 2 ; //length of the output vector
    NumericVector dists(nout); //preallocating output vector
    for (int i= 0 ,k= 0 ;i<nin- 1 ;i++) {
        double a = v[i];
        for (int j=i+ 1 ;j<nin;j++,k++) {
            double b = v[j];
            dists[k] = abs(a-b);
        }
    }
    return dists;
}

实验

我用10.000个随机数的向量分别运行了3次不同的实现（避免冗余时需要49.995。000距离计算，否则要使用100.000.000），并计算了CPU时间和内存消耗的中位数。 我使用了2017 MacBook Pro，2.3 GHz Intel Core i5、16 GB RAM，运行带有R 3.6.0，Julia 1.0.3和XCode 10.1的Mac OS 10.14。

正如预期的那样，R中基于循环的实现是最慢的（并且在提供JIT的3.4版之前会慢得多）。 通过矢量化，我们减少了计算时间，但增加了内存消耗，随着输入大小的增加，这可能成为问题。 即使进行了向量化工作，我们仍然离R的dist函数的性能还差得远。

Rcpp允许减少计算时间和内存需求，胜过R的核心实现。 这并不奇怪，因为R的dist函数更加灵活，添加了多个选项和输入验证。

可以将C / CPP代码注入R脚本虽然很棒，但是现在我们正在使用两种编程语言，并且已经失去了CPP代码交互式编程的好处。

Julia的方式

Julia的编程思想与R完全不同。 最有效的解决方案是通过一种基于循环的方法来预分配内存。

function distances_among_pairs_loop(v)
    nin = length(v)
    nout = div(nin * nin - nin, 2 ) #length of the output vector
    dists = Vector {eltype(v)}(undef,nout) #preallocating output
    k = 1 #current position in the output vector
    for i in 1 :(nin- 1 )
        a = v[i]
        for j in (i+ 1 ):nin
            b = v[j]
            dists[k] = abs(a-b)
            k += 1
        end
    end
    dists
end

如果您想编写更少的代码，则可以以降低计算效率为代价。 压缩代码的一种方法是通过理解。

distances_among_pairs_comp(v) = 
    [abs(v[i]-v[j])for i in eachindex(v), j in eachindex(v) if i<j]

上面的理解比基于循环的实现更紧凑，同时体现了相同的逻辑。 这种方法的主要缺点是输出矢量未预先分配。 由于编译器无法预测输出的大小，因此输出会根据需要动态增长。

最后，如果您喜欢向量化方法，那么Julia中也可以选择这种方法。 基于外部函数转换R的实现如下所示：

function distances_among_pairs_outer(v)
    dists = v .- v'
    abs.(dists[tril!(trues(size(dists)), - 1 )])
end

与R中不同，我们希望这是效率最低的方法，因为它需要更多的内存和不必要的（冗余）操作。

结果

Julia通过开箱即用地提供类似C的性能而脱颖而出。 代码紧凑性和效率之间的权衡非常清楚，类似C的代码可提供类似C的性能。 理解是一个很好的折衷方案，因为它们更易于编写代码，不易出现错误，并且在此问题上的内存效率相同。

当使用R解决计算密集型任务时，我们希望找到一个专门的函数来解决我们的问题。 如果没有专用功能，我们要么需要在C / CPP中对其进行编程，要么需要通过向量化途径进行。 如果选择第二个，则最终结果可能会与第一个结果相差甚远。

总结思想

向量化，理解，地图过滤器减少都是非常棒的工具，可以节省您的时间并提供更紧凑，更易读的代码。 但是，由于编程语言的性能限制，不应强迫程序员使用它们。

如果有直接而有效的实施方式，您不想花费时间去尝试几种方法来实施解决方案。

Julia允许您在代码紧凑性和计算效率之间进行选择，从而使这种权衡非常明确。 您可以实现类似C的循环，类似R的矢量化或类似Python的理解。

由你决定。 您可以优化编写解决方案所需的时间与解决方案运行所花费的时间之间的比率。 您可以实现自己的算法，而无需依赖第二种语言。

Julia（Julia）比Python和R年轻得多，他正在通过数据科学界奋斗。 13年前我开始使用Matlab时就非常喜欢，因为我可以与数据进行交互。 我现在喜欢如何更快地编写更快的代码，以及如何自由地在Julia中实现几乎任何算法。 放手，我相信您也会喜​​欢它！

代码可在 github.com/dcmoura/blogposts获得

原始出版物@ 迈向数据科学

您可以在Tweeter ， LinkedIn和Vimeo上找到我（请查看我的datavizs！）。

翻译自: https://hackernoon.com/freeing-the-data-scientist-mind-from-the-curse-of-vectorization-julia-to-the-rescue-0c3z308v

栅格数据和矢量数据矢量化