Julia中文社区

common-random-variables

动机 当我们模拟带随机成分的模型时,它通常是有利于分解成如下面这样的结构:

  1. 含有需要我们来描述结构关系的参数θ,
  2. 含有独立于参数的噪音 ε,
  3. 含有能够生成可观察数据或矩阵的函数 f(θ,ε)

这样做可以让我们使用常见的随机变量(也就是常见的随机数),这是一种常见技巧能简单地固定ε,尝试不同的θ。这样做可使函数f成为可微分的, 使f能用基于导数的优化算法(例如,最大似然法或MAP)或基于MCMC方法的衍生方法。并且它还用于减少模拟矩的方差。

在Julia中实现这种技术时,我经常创建一个能保存参数ε的盒子结构,它会分配一个需原地更新的数组空间。由于 DynamicHMC.jl 的 再版设计在出现时会以更加规范的方式容纳模拟似然方法,所以我想探索其他选择,其中最重要的是StaticArrays.jl

在这里,我会用一个简单的玩具模型对Julia v0.6.2的替代品进行基准测试。对比TL,DR来说:StaticArrays.jl 快150倍, 并且这里不依赖于所用的静态数组的可变性,甚至每次重新分配 时都会产生新的ε 。

设置问题: 简单设置一下:我们准备出$M$个观察对象,然后噪音是:

εi,j∼N(0,1)for i=1,…,M;j=1,…,7. 我们的参数是μ和σ,对于每个i我们来算

这是非线性变换后的样本平均值。然后这个7是偶然选取的,它来自简化我正在处理的实际问题。我们对Ai的样本平均值感兴趣。我刻意避免只是微优化每个版本,这样才能反映出怎么处理现实生活中的问题。 我们将这个通用接口编代码如下: using BenchmarkTools using StaticArrays using Parameters

common interface

"Dimension of noise ϵ for each observation." const EDIM = 7

«»" Common random variables. The user needs to define

  1. observation_moments, which should use observation_moment,
  2. newcrv = update!(crv), which returns new common random variables, potentially (but not necessarily) overwriting crv. """ abstract type CRVContainer end

observation_moment(ϵ, μ, σ) = mean(@. exp(μ + σ * ϵ))

average_moment(crv::CRVContainer, μ, σ) = mean(observation_moments(crv, μ, σ))

"Calculate statistics, making N draws, updating every Lth time." function stats(crv, μ, σ, N, L) _stat() = (N % L == 0 && (crv = update!(crv)); average_moment(crv, μ, σ)) [_stat() for _ in 1:N] end

我写stats函数的方式代表了我如何使用HMC / NUTS:在相同轨迹上的模拟力矩就用相同的ε来计算,ε也会根据每一条轨迹来更新。当然,参数会沿轨迹变化,不过这里它们没有变化,但这并不影响基准。 update! 函数的语义能够直接修改功能和风格。

使用一个预分配矩阵 这是我写代码中用到的主要方法。1 ε在矩阵中的列里,将它们映射为片状更好,然后在用到 observation_moment时会被映射。 update! 函数重写了这些内容。 ""« Common random variables are stored in columns of a matrix. »"" struct PreallocatedMatrix{T} <: CRVContainer ϵ::Matrix{T} end

PreallocatedMatrix(M::Int) = PreallocatedMatrix(randn(EDIM, M))

update!(p::PreallocatedMatrix) = (randn!(p.ϵ); p)

observation_moments(p::PreallocatedMatrix, μ, σ) = vec(mapslices(ϵ -> observation_moment(ϵ, μ, σ), p.ϵ, 1)) 利用静态向量 我们在各种静态向量实现间共享了以下内容: "Common random variables as a vector of vectors, in the ϵs." abstract type CRVVectors <: CRVContainer end

observation_moments(p::CRVVectors, μ, σ) = map(ϵ -> observation_moment(ϵ, μ, σ), p.ϵs) 我发现在这里使用映射比使用上面说的映射切片要更加直观。 静态向量,预分配容器 在使用静态向量的第一个版本中, 我们将 SVector 保留在 Vector中,并适时更新。. struct PreallocatedStaticCRV{K, T} <: CRVVectors ϵs::Vector{SVector{K, T}} end

PreallocatedStaticCRV(M::Int) = PreallocatedStaticCRV([@SVector(randn(EDIM)) for _ in 1:M])

function update!(p::PreallocatedStaticCRV) @unpack ϵs = p @inbounds for i in indices(ϵs, 1) ϵs[i] = @SVector(randn(EDIM)) end p end 可变静态向量与重写 我们修改这个来使用可变向量,这应该没有什么区别。 struct MutableStaticCRV{K, T} <: CRVVectors ϵs::Vector{MVector{K, T}} end

MutableStaticCRV(M::Int) = MutableStaticCRV([@MVector(randn(EDIM)) for _ in 1:M])

function update!(p::MutableStaticCRV) @unpack ϵs = p @inbounds for i in indices(ϵs, 1) randn!(ϵs[i]) end p end 每次分配的静态向量 最后,这就是第三版, struct GeneratedStaticCRV{K, T} <: CRVVectors ϵs::Vector{SVector{K, T}} end

GeneratedStaticCRV(M::Int) = GeneratedStaticCRV([@SVector(randn(EDIM)) for _ in 1:M])

update!(p::GeneratedStaticCRV{K, T}) where {K, T} = GeneratedStaticCRV([@SVector(randn(T, K)) for _ in indices(p.ϵs, 1)]) 结果 运行如下代码: @btime mean(stats(PreallocatedMatrix(100), 1.0, 0.1, 100, 10)) @btime mean(stats(PreallocatedStaticCRV(100), 1.0, 0.1, 100, 10)) @btime mean(stats(MutableStaticCRV(100), 1.0, 0.1, 100, 10)) @btime mean(stats(GeneratedStaticCRV(100), 1.0, 0.1, 100, 10)) 我们获得:

作为将来改进的预览,我在当前的master 上尝试了 PreallocatedMatrix  (master 将成为Julia v0.7),这获得了3.5 ms (2.46 MiB)的结果, 这正是我们希望看到的。2 结论是, StaticArrays 简化了逻辑,同时加速了代码的速度.。我特别喜欢上一版本(GeneratedStaticCRV), 因为它不需要事先考虑类型。. 虽然这里的例子很简单,但实际上我会使用自动微分,这使得提前确定缓冲区类型更具挑战性。. 我希望将来我会转换到更“无缓冲”的风格,并相应地为DynamicHMC.jl设计界面。

我们为什么创造Julia语言

这是一篇重译,参考了2012年的一篇豆瓣(链接在最后)。

简短来讲,是因为我们很贪婪。

我们之中有些是使用MATLAB的重量级用户,有些是来自Lisp的极客,还有一些是来自Python和Ruby的魔法师,甚至还有来自Perl社区的大魔法师。我们之中还有从胡子都没长齐时就开始使用Mathematica的。其中的有些人现在都没长胡子喱!我们像是疯了一样用R画了越来越多的图,而C是我们的硬核摇滚(也有大杀器之意)。

我们热爱所有这些语言,他们实在很好很强大。在我们从事的领域(科学计算,机器学习,数据挖掘,大规模线性代数计算,分布式和并行计算)中,每一种语言都对某一项工作的一项特定需求非常完美,但是却无法胜任其它需求。于是使用什么语言都需要我们去权衡。

而我们很贪婪,我们还想要更多。

我们想要的是一个自由开源的语言,并且它同时拥有C的速度和Ruby的动态性;我们想要一个具有同像性(可以将语言的脚本本身当作数据进行处理)的语言,它有着真正的和lisp一样的宏,但是却像Matlab一样有着显然的,类似于数学表达式的标记;我们想要一个既可以像Python一样作为通用编程语言的工具,又可以像R那样适用于统计分析,能像Perl那样自然地处理字符串,能像Matlab那样给力地处理矩阵运算,它还要能像shell一样作为胶水将各种程序粘合在一起;我们想要一个简单易学的语言,同时它还能让最苛刻的魔法师们(hackers)开心。我们希望它是交互式的,但我们也希望它能被编译。

(我们刚刚有提它要和C一样快嘛?!)

当我们在构思这些需求的时候,我们发现它还得有Hadoop这样强大的分布式能力,却不想要Hadoop里面那些冗长Jave和XML模板,更不想被被迫在几个GB的日志文件和几百台机器里找bug。我们不想要那些令人费解的层次结构。我们想让简单的标量循环能被编译成仅用寄存器和一块CPU的干净的机器码。我们希望简单地写下A*B就能够在成千上万的机器上用成千上万地运算来计算这个庞大的矩阵乘法。

如无必要,那就不用声明类型。但当我们需要多态函数(polymorphic functions)时,我们也想要用泛型编程(generic programming)仅仅书写一次算法,就能够在无限多的类型上使用。我们想要多重派发(multiple dispatch)来为一个函数所有可能的参数选出最佳的执行方法。这些参数可能有着不同定义,不同类型,但是却有着相同功能。在拥有以上能力的同时,我们还希望这种语言简单,干净。

要求有点多,是不是?

尽管我们意识到了自己有多贪心,我们还是想要拥有这些功能。大概在两年半之前,我们开始创造这种能满足我们贪念的语言。它还没有完工——但是已经可以发布一个1.0版本了(其实等了6年才要发布)——我们创造的这个语言叫做Julia。它已经实现了我们这次乱七八糟需求的90%,而现在她需要来自更多人的乱七八糟的需求,来让她走得更远。如果你也是一位贪心不足,不可理喻,需索无度的码场二逼青年,希望你能来试试这个东东。

作者:Jeff Bezanson, Stefan Karpinski, Viral Shah, Alan Edelman
译者:Roger
翻译自:Why we create Julia
参考自:豆瓣:为什么我们要创造Julia