Julia程序效率介绍

编辑:李东风

2020-04-05

Julia的语言特点

Julia语言是一种历史很短的计算机语言, 公开发布于2012年。 其设计理念就是希翼兼有Python、R、Matlab这样的动态语言的易用性, 以及C、C++、Java这样的静态语言的运行速度。 所以Julia很适合用来做统计和金融计算。

Julia语言的特点有:

  • 动态语言;
  • 使用基于LLVM的动态编译技术,可以动态生成高效的运行代码;
  • 不需要声明类型也一般能够产生高效代码, 需要时声明类型可以提高效率。

程序性能评测

程序计时

@time宏对一个表达式计时和查看内存分配情况,如:

In [3]:
@time sqrt.(rand(10_000_000));
  0.118549 seconds (4 allocations: 152.588 MiB, 4.92% gc time)

结果中有程序运行时间的秒数, 进行了10次内存分配, 总计分配的内存数(中间可能也有释放的内存)。

@time宏实行一个表达式, 显示运行时间、分配内存、垃圾收集时间比例, 然后返回表达式的结果。

@time也可以对一个结构计时,如

In [4]:
s = 0.0
@time for i=1:10_000_000
    global s
    s += sqrt(rand())
end
  1.236361 seconds (20.00 M allocations: 305.176 MiB, 13.81% gc time)

这个程序慢了很多,分配内存的次数和总量也比较多。

还有一个@timev宏与@time类似, 提供的信息更多。 @time@timev仅是显示表达式实行的时间和内存分配信息, 不影响表达式返回结果, 所以在程序中使用这两个函数计时一般不影响程序的结果。

程序瓶颈查找(profiler)

为了提高程序效率, 并不是所有的程序代码都需要仔细地考察并改进效率。 大部分的程序可能仅仅实行一两次, 可能仅需要实行几秒钟到几分钟。 这样的程序不需要想办法提高实行效率, 只要程序结果正确就可以了。

有些程序需要实行许多次, 可能会提供给许多用户使用, 这样的程序即使耗时并不多, 也有必要去设法改进其运行效率。 有些程序虽然不需要实行许多次, 但是运行所需的时间可能很长, 比如几小时到几天, 这样的程序需要设法改进效率。

对一段比较长的程序, 如果需要改进效率, 并不需要考虑所有的代码行, 而是只要改进程序实行消耗时间最多的部分就可以了, 这样的部分称为程序的瓶颈。 Julia内置了一个profiler工具, 用定时抽查的方式查看那些代码行被实行得最多。 这样的工具在使用时不需要修改源程序, 对运行速度的影响也很小。

因为Julia使用了即时编译, 所以一个函数在第一次被调用时需要有额外的编译时间, 第二次调用就没有这个问题了。 在用profiler分析一个函数是应该从第二次调用开始分析。

在实行程序之前,调用:

using Profile
Profile.clear()

然后实行要分析的程序。 实行完后调用:

Profile.print()

将显示程序的实行次数, 包括调用的库函数的实行次数, 实行次数多的程序语句是需要优化的。

ProfileView包可以提供一个交互的图形显示, 直观地显示哪些语句耗时多。 编辑个人认为这些工具现在显示的信息还过于繁琐, 不利于一般用户使用。

BenchmarkTools

为了对运行时间较短的程序计时, 通常会人为地重复该程序。 BenchmarkTools提供了更完善的计时统计, 以及系统的性能比较(benchmark)能力。

例如,如下的简短程序:

In [5]:
A = rand(2); B = rand(2); C = zeros(2);
C = A + B
Out[5]:
2-element Array{Float64,1}:
 1.254767738514001
 1.6831627540750678

对第二行计时如下:

In [6]:
using BenchmarkTools
@benchmark C = A + B
Out[6]:
BenchmarkTools.Trial: 
  memory estimate:  96 bytes
  allocs estimate:  1
  --------------
  minimum time:     91.508 ns (0.00% GC)
  median time:      99.056 ns (0.00% GC)
  mean time:        156.384 ns (2.54% GC)
  maximum time:     2.817 μs (88.40% GC)
  --------------
  samples:          10000
  evals/sample:     954

换用加点写法:

In [7]:
@benchmark C .= A .+ B
Out[7]:
BenchmarkTools.Trial: 
  memory estimate:  48 bytes
  allocs estimate:  2
  --------------
  minimum time:     375.995 ns (0.00% GC)
  median time:      386.505 ns (0.00% GC)
  mean time:        411.862 ns (0.57% GC)
  maximum time:     12.348 μs (96.51% GC)
  --------------
  samples:          10000
  evals/sample:     200

这个版本性能更差。

换用显式循环:

In [8]:
@benchmark for i=1:2
    C[i] = A[i] + B[i]
end
Out[8]:
BenchmarkTools.Trial: 
  memory estimate:  96 bytes
  allocs estimate:  6
  --------------
  minimum time:     246.130 ns (0.00% GC)
  median time:      272.755 ns (0.00% GC)
  mean time:        298.867 ns (0.48% GC)
  maximum time:     2.870 μs (90.62% GC)
  --------------
  samples:          10000
  evals/sample:     323

不如向量化版本, 但优于加点版本。

以上的程序是在全局名字空间(命令行对应的名字空间)中实行的。 如果调用函数, 结果则不同:

In [9]:
A = rand(2); B = rand(2); C = zeros(2);

function vecadd1(a,b,c)
    c = a + b
    return nothing
end

function vecadd2(a,b,c)
    c .= a .+ b
    return nothing
end

function vecadd3(a,b,c)
    for i=1:length(c)
        c[i] = a[i] + b[i]
    end
    return nothing
end
@benchmark vecadd1(A, B, C)
Out[9]:
BenchmarkTools.Trial: 
  memory estimate:  96 bytes
  allocs estimate:  1
  --------------
  minimum time:     81.744 ns (0.00% GC)
  median time:      83.091 ns (0.00% GC)
  mean time:        91.030 ns (3.67% GC)
  maximum time:     1.884 μs (95.15% GC)
  --------------
  samples:          10000
  evals/sample:     964
In [10]:
@benchmark vecadd2(A, B, C)
Out[10]:
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     39.717 ns (0.00% GC)
  median time:      40.423 ns (0.00% GC)
  mean time:        41.850 ns (0.00% GC)
  maximum time:     127.117 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     992
In [11]:
@benchmark vecadd3(A, B, C)
Out[11]:
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     33.802 ns (0.00% GC)
  median time:      34.003 ns (0.00% GC)
  mean time:        35.109 ns (0.00% GC)
  maximum time:     106.338 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     994

加点和显式循环的版本性能相近, 都优于向量化版本, 比命令行版本也提高了很多。 这可能与即时编译系统对函数的优化有关。

数据类型与运行效率

类型声明

Julia允许对函数自变量、函数返回值、函数的局部变量声明数据类型, 但是一般不需要声明也能产生运行效率很好的实行代码。 声明方式是在变量名后面用双冒号::后面跟类型名,如

In [12]:
function f_tp01(x::Float64) :: Float64
    local y::Float64
    y = (1 + x)^2
    y = sin(y)
    return y
end
f_tp01(1.5)
Out[12]:
-0.03317921654755682

说明类型有可能会提高程序效率, 但是也限制了函数调用方法, 尤其是整型与实型会看成两种不同类型,所以:

f_tp01(1)

结果:

MethodError: no method matching f_tp01(::Int64)
Closest candidates are:
  f_tp01(!Matched::Float64) at In[12]:2

Stacktrace:
 [1] top-level scope at In[13]:1

上面的调用错误是因为f_tp01()函数只允许浮点型自变量, 而1不是浮点型, Julia在调用函数时不会自动将整型转换成浮点型。 因为Julia的四则运算的设计是允许将整型与浮点型混合运算的, 所以四则运算不存在这样的问题。

函数定义时如果不声明自变量类型, 则不存在这个问题,如:

In [14]:
function f_tp01b(x)
    y = (1 + x)^2
    y = sin(y)
    return y
end
f_tp01b(1.5)
Out[14]:
-0.03317921654755682
In [15]:
f_tp01b(1)
Out[15]:
-0.7568024953079282

所以自定义函数中如非必要并不是所有自变量、返回值、局部变量都声明类型就更好。

为了使得f01()能够对其它数也适用, 可以写一个更一般的类型作为自变量然后转换成浮点型,如:

In [16]:
function f_tp01(x::Number) :: Float64
    local y::Float64
    y = (1 + Float64(x))^2
    y = sin(y)
    return y
end
f_tp01(1)
Out[16]:
-0.7568024953079282

多重派发

Julia的同一个函数可以根据自变量个数和自变量类型的不同而实行不同的操作, 称为“多重派发”(multiple dispatch)。 这有些像是面向对象语言的方法重载(overloading), 但是:

  • 方法重载是编译时根据对象的类决定的;
  • 方法重载仅能根据其依附的对象而选择对应的方法,即所谓单重派发;
  • Julia是动态语言,其函数调用是运行时动态决定的;
  • Julia可以根据其多个自变量的类型而选取不同方法。

多重派发可以让多个类似的操作共享同一个函数名。 比如, 可以这样设计绘图函数plot(), 可以仅输入一个数值型向量y, 作序列图; 可以输入两个数值型向量xy,作散点图; 可以输入一个字符串向量x和数值型向量y, 作条形图。 不同的自变量个数和类型决定了要实行的操作, 这些操作称为plot()函数的“方法”。

类型稳定性

Julia不要求声明函数自变量、返回值和局部变量类型, 但是如果函数返回值不能被其输入数据类型确定, 这样的函数很难得到高效的实行代码。 这种要求称为程序的“类型稳定性”。

例如,下面的程序:

In [18]:
function f_tp02(x)
    if x >= 0
        return sqrt(x)
    else
        return "Negative number have no square root!"
    end
end
Out[18]:
f_tp02 (generic function with 1 method)

这个函数从表面上看很正常, 但是违反了类型稳定性。 因为不论输入x为什么数, 返回值的类型都依赖于x的数值而不能由x的数据类型决定, 这就使得依赖于返回值的代码无法确定返回值类型从而无法优化。

用宏@code_warntype可以检查一个函数调用中的类型稳定性问题,如

In [19]:
@code_warntype f_tp02(1.5)
Variables
  #self#::Core.Compiler.Const(f_tp02, false)
  x::Float64

Body::Union{Float64, String}
1 ─ %1 = (x >= 0)::Bool
└──      goto #3 if not %1
2 ─ %3 = Main.sqrt(x)::Float64
└──      return %3
3 ─      return "Negative number have no square root!"

结果中的类型不稳定问题会用红色显示。

另外,如果一个函数反复被调用或者内部有很多次循环, 某些局部变量在多次调用或者多次循环中不能预先确定类型或者类型被改变, 也会造成程序性能急剧下降。 如:

In [20]:
function f_tp03(x)
  s = 0
  for xi in x
    s += sqrt(xi)
  end 
  return s
end
Out[20]:
f_tp03 (generic function with 1 method)

变量s开始是整型,但在循环中变成了浮点型。初始化可写成s = 0.0

函数的效率

由于Julia的多重派发、动态编译特性, Julia程序最好分解成短小、输入类型明确的多个函数。 函数调用的成本很低, 动态编译可以将函数编译成针对特定输入的高效代码。

全局变量的问题

为了运行效率考虑, 应尽量避免使用全局变量。 全局变量的取值乃至于类型都可能被不同位置的代码改变, 使得编译器很难预估全局变量的类型, 造成编译代码低效。 尽可能使用局部变量与函数自变量保存和传递数据。 Julia的函数是按引用进行参数传递的, 所以数组的传递不会制作副本, 没有内存和时间的额外开销。

Julia的编译器优化主要是针对函数, 所以对效率敏感的代码都应该写成函数而不是直接在全局变量空间实行。

const声明的全局变量不能改变类型, 这样的全局变量对程序效率的影响较小。

考虑下面的利用了全局变量作为选项的程序:

In [21]:
POWER=2
function f_gl01(x::Vector{Float64})
    s = 0.0
    for xi in x
        s += xi ^ POWER
    end
    return s
end
f_gl01(rand(10))

using BenchmarkTools
@benchmark f_gl01(rand(100_000))
Out[21]:
BenchmarkTools.Trial: 
  memory estimate:  5.34 MiB
  allocs estimate:  300002
  --------------
  minimum time:     7.114 ms (0.00% GC)
  median time:      8.045 ms (0.00% GC)
  mean time:        8.918 ms (1.96% GC)
  maximum time:     23.844 ms (0.00% GC)
  --------------
  samples:          560
  evals/sample:     1

将上面的POWER声明为常数可以提高效率。 Julia用const关键字声明全局变量常数, 但是这种全局变量仅仅是不允许改变数据类型, 变量值还是允许改变的, 改变变量值会有警告信息。

修改后的程序如:

In [22]:
const P_const=2
function f_gl01b(x::Vector{Float64})
    s = 0.0
    for xi in x
        s += xi ^ P_const
    end
    return s
end
f_gl01b(rand(10))

using BenchmarkTools
@benchmark f_gl01b(rand(100_000))
Out[22]:
BenchmarkTools.Trial: 
  memory estimate:  781.33 KiB
  allocs estimate:  2
  --------------
  minimum time:     266.599 μs (0.00% GC)
  median time:      614.299 μs (0.00% GC)
  mean time:        700.078 μs (10.89% GC)
  maximum time:     7.292 ms (91.53% GC)
  --------------
  samples:          7108
  evals/sample:     1

比原来效率提高了一个数量级。

在调用全局变量的位置标识其类型也可以提高效率,如

In [23]:
POWER=2
function f_gl01c(x::Vector{Float64})
    s = 0.0
    for xi in x
        s += xi ^ POWER::Int
    end
    return s
end
f_gl01c(rand(10))

using BenchmarkTools
@benchmark f_gl01c(rand(100_000))
Out[23]:
BenchmarkTools.Trial: 
  memory estimate:  781.33 KiB
  allocs estimate:  2
  --------------
  minimum time:     907.800 μs (0.00% GC)
  median time:      1.229 ms (0.00% GC)
  mean time:        1.291 ms (4.95% GC)
  maximum time:     6.446 ms (68.08% GC)
  --------------
  samples:          3862
  evals/sample:     1

当然,这个问题实际上完全不需要使用全局变量, 只要作为函数自变量即可:

In [24]:
function f_gl01d(x::Vector{Float64}, p)
    s = 0.0
    for xi in x
        s += xi ^ p
    end
    return s
end
f_gl01d(rand(10), 2)

using BenchmarkTools
@benchmark f_gl01d(rand(100_000), 2)
Out[24]:
BenchmarkTools.Trial: 
  memory estimate:  781.33 KiB
  allocs estimate:  2
  --------------
  minimum time:     265.100 μs (0.00% GC)
  median time:      599.699 μs (0.00% GC)
  mean time:        631.232 μs (10.43% GC)
  maximum time:     7.973 ms (90.86% GC)
  --------------
  samples:          7888
  evals/sample:     1

@code_warntype分析f_gl01()的问题:

In [25]:
@code_warntype f_gl01(rand(10))
Variables
  #self#::Core.Compiler.Const(f_gl01, false)
  x::Array{Float64,1}
  s::Any
  @_4::Union{Nothing, Tuple{Float64,Int64}}
  xi::Float64

Body::Any
1 ─       (s = 0.0)
 %2  = x::Array{Float64,1}
       (@_4 = Base.iterate(%2))
 %4  = (@_4 === nothing)::Bool
 %5  = Base.not_int(%4)::Bool
└──       goto #4 if not %5
2 ┄ %7  = @_4::Tuple{Float64,Int64}::Tuple{Float64,Int64}
       (xi = Core.getfield(%7, 1))
 %9  = Core.getfield(%7, 2)::Int64
 %10 = s::Any
 %11 = (xi ^ Main.POWER)::Any
       (s = %10 + %11)
       (@_4 = Base.iterate(%2, %9))
 %14 = (@_4 === nothing)::Bool
 %15 = Base.not_int(%14)::Bool
└──       goto #4 if not %15
3 ─       goto #2
4 ┄       return s

输出中的红色Any都提示有类型不确定问题。修改然后再用@code_warntype分析:

In [26]:
POWER=2
function f_gl01e(x::Vector{Float64})::Float64
    s::Float64 = 0.0
    for xi in x
        s += xi ^ POWER::Int
    end
    return s
end
f_gl01e(rand(10))

using BenchmarkTools
@benchmark f_gl01e(rand(100_000))
Out[26]:
BenchmarkTools.Trial: 
  memory estimate:  781.33 KiB
  allocs estimate:  2
  --------------
  minimum time:     873.500 μs (0.00% GC)
  median time:      1.243 ms (0.00% GC)
  mean time:        1.317 ms (4.86% GC)
  maximum time:     9.083 ms (77.60% GC)
  --------------
  samples:          3786
  evals/sample:     1
In [26]:
@code_warntype f_gl01e(rand(10))
Variables
  #self#::Core.Compiler.Const(f_gl01e, false)
  x::Array{Float64,1}
  s::Float64
  @_4::Union{Nothing, Tuple{Float64,Int64}}
  xi::Float64

Body::Float64
1 ─ %1  = Main.Float64::Core.Compiler.Const(Float64, false)
 %2  = Base.convert(Main.Float64, 0.0)::Core.Compiler.Const(0.0, false)
       (s = Core.typeassert(%2, Main.Float64))
 %4  = x::Array{Float64,1}
       (@_4 = Base.iterate(%4))
 %6  = (@_4 === nothing)::Bool
 %7  = Base.not_int(%6)::Bool
└──       goto #4 if not %7
2 ┄ %9  = @_4::Tuple{Float64,Int64}::Tuple{Float64,Int64}
       (xi = Core.getfield(%9, 1))
 %11 = Core.getfield(%9, 2)::Int64
 %12 = s::Float64
 %13 = xi::Float64
 %14 = Core.typeassert(Main.POWER, Main.Int)::Int64
 %15 = (%13 ^ %14)::Float64
 %16 = (%12 + %15)::Float64
 %17 = Base.convert(Main.Float64, %16)::Float64
       (s = Core.typeassert(%17, Main.Float64))
       (@_4 = Base.iterate(%4, %11))
 %20 = (@_4 === nothing)::Bool
 %21 = Base.not_int(%20)::Bool
└──       goto #4 if not %21
3 ─       goto #2
4 ┄ %24 = Base.convert(%1, s)::Float64
 %25 = Core.typeassert(%24, %1)::Float64
└──       return %25

结果中没有需要改进的代码了。

行内函数(inline functions)

Julia编译器自动将较小的类型明确的函数调用变成行内函数, 即被调用的函数的实际代码被转移到调用者内部。 行内函数提高了程序效率但是会使得程序变大。

Julia自动确定哪些函数变成行内函数, 然而, 有些反复调用的代码, 尤其是用在循环和内存循环的代码, 比如数组用下标访问元素的函数, 最好人为地规定成行内函数。 用@inline宏说明一个函数定义就可以规定其变成行内函数。 如

@inline function f(x)
    ...
end

用宏提高效率

宏(macro)是类似C语言预处理的语言构件。 宏可以生成程序代码, 可以对程序代码进行操作, 比如前面的@time@code_warntype等。

宏可以用来对重复性的代码进行简写, 也可以通过将运行时才实行的任务提前到代码编译阶段完成从而提高效率。

以多项式求值的函数为例。 多项式的系数按升序保存在一维数组中。 按照通常的写法,程序为:

In [27]:
function f_ma01(x, a...)
    p = zero(x)
    for i = 1: length(a)
        p += a[i] * x ^(i-1)
    end
    return p
end
Out[27]:
f_ma01 (generic function with 1 method)

这里没有用秦九韶法进行算法优化,直接用了幂函数。 对多项式$f(x) = 1 + 2 x + 3 x^2 + 4 x^3$, 可以用f_ma01()定义如下的多项式函数:

In [28]:
f_ma01b(x) = f_ma01(x, 1, 2, 3, 4)
Out[28]:
f_ma01b (generic function with 1 method)

调用如

In [29]:
f_ma01b(1.0)
Out[29]:
10.0
In [30]:
f_ma01b(1.5)
Out[30]:
24.25
In [31]:
f_ma01b(1)
Out[31]:
10

测试程序效率如下:

In [32]:
@benchmark f_ma01b(1.5)
Out[32]:
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     1.899 ns (0.00% GC)
  median time:      2.000 ns (0.00% GC)
  mean time:        2.009 ns (0.00% GC)
  maximum time:     13.700 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1000

将程序改成用秦九韶法, 比如$f(x) = 1 + 2 x + 3 x^2 + 4 x^3$ 写成 $$ f(x) = 1 + x( 2 + x (3 + 4 x)) $$ 这样可以减少乘法运算的次数。 算法每一步将前一步的结果乘以x, 加上降序的后一个系数。 将f_ma01()用这种方法改写:

In [33]:
function f_ma02(x, a...)
    p = zero(x)
    for i = length(a):-1:1
        p = a[i] + p*x
    end
    return p
end
Out[33]:
f_ma02 (generic function with 1 method)

多项式$f(x) = 1 + 2 x + 3 x^2 + 4 x^3$ 的函数为:

In [34]:
f_ma02b(x) = f_ma02(x, 1, 2, 3, 4)
f_ma02b(1.5)
@benchmark f_ma02b(1.5)
Out[34]:
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     1.899 ns (0.00% GC)
  median time:      2.000 ns (0.00% GC)
  mean time:        1.980 ns (0.00% GC)
  maximum time:     34.300 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1000

两个函数性能相同, 可能是Julia本身进行了优化。

这个算法还是不够优化, 因为每次调用f_ma02()时都需要访问系数数组, 而实际的多项式的系数是常数, 所以, 直接按$f(x) = 1 + x( 2 + x (3 + 4 x))$计算效率会更高,如:

In [35]:
f_ma03(x) = 1 + x*(2 + x*(3 + 4 * x))
f_ma03(1.5)
@benchmark f_ma03(1.5)
Out[35]:
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     0.001 ns (0.00% GC)
  median time:      0.001 ns (0.00% GC)
  mean time:        0.034 ns (0.00% GC)
  maximum time:     0.101 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1000

第三个版本效率提升了三个数量级。 问题是f_ma03()函数是完全根据多项式阶数与系数手写代码实现的, 大家需要就这个手写代码的过程用Julia程序完成, 这就是宏的作用。

一般的生成一个多项式函数的宏如下:

In [36]:
macro mac_01(x, a...)
    ex = esc(a[end])
    for i = length(a)-1:-1:1
        ex = :(muladd(t, $ex, $(esc(a[i]))))
    end
    Expr(:block, :(t = $(esc(x))), ex)
end
Out[36]:
@mac_01 (macro with 1 method)

生成多项式函数如:

In [37]:
f_ma04(x) = @mac_01(x, 1, 2, 3, 4)
Out[37]:
f_ma04 (generic function with 1 method)
In [38]:
f_ma04(1.5)
@benchmark f_ma04(1.5)
Out[38]:
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     0.001 ns (0.00% GC)
  median time:      0.001 ns (0.00% GC)
  mean time:        0.036 ns (0.00% GC)
  maximum time:     1.800 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1000

程序效率与手写的优化代码相当。

派生函数(generated functions)是与宏类似的函数定义方式, 可以起到类似手写优化代码的效果, 可以将原来一般性的循环写成直接的表达式。

关键字参数问题

Julia的多重派发是仅根据位置参数, 包括没有缺省值的和有缺省值的位置参数的类型来进行的; 关键字参数不影响多重派发, 这样, 用到关键字参数的函数就无法得到类型信息的好处, 效率相对而言比较低。 关键字参数适用于有很多选项的用户直接调用的参数, 如果是反复循环的对性能要求高的函数则应避免使用关键字参数。 例如:

In [39]:
f_kw01(x; y=2, z=2) = x^y + x^z
f_kw01(1.5, y=2, z=3)
@benchmark f_kw01(1.5, y=2, z=3)
Out[39]:
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     2.299 ns (0.00% GC)
  median time:      2.399 ns (0.00% GC)
  mean time:        2.368 ns (0.00% GC)
  maximum time:     34.901 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1000
In [40]:
f_kw01b(x, y=2, z=2) = x^y + x^z
f_kw01b(1.5, 2, 3)
@benchmark f_kw01b(1.5, 2, 3)
Out[40]:
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     2.599 ns (0.00% GC)
  median time:      2.699 ns (0.00% GC)
  mean time:        2.706 ns (0.00% GC)
  maximum time:     35.401 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1000

两个版本性能相近, 可能是Julia做了优化。

基本计算

整数的表示与溢出

Julia的数值直接使用硬件的数值表示。 整数值的计算不进行溢出检查, 所以利用整数值计算时需要用户自己注意溢出问题,如:

In [41]:
2^62
Out[41]:
4611686018427387904
In [42]:
2^63
Out[42]:
-9223372036854775808

不进行溢出检查是为了程序效率必须付出的代价。

如果确实有可能发生溢出, Int128可以保存更多位, 进一步可以改用浮点数Float64, 或者用无限位数的整数类型BigInt, BigInt效率会很差。 如

In [43]:
BigInt(2)^63
Out[43]:
9223372036854775808

即使用BigInt()说明, 也仅是将输入转换成了无限精度整型, 转换之前的溢出是没有控制的,比如:

In [44]:
BigInt(2^63)
Out[44]:
-9223372036854775808

要注意, 直接写的整数根据数字位数, 较短时为Int类型, 超过Int范围首先选用Int64和Int128, 超过Int128后自动认为是BigInt类型。 如:

In [45]:
typeof(123456789012345678)
Out[45]:
Int64
In [46]:
typeof(12345678901234567890)
Out[46]:
Int128
In [47]:
typeof(1234567890123456789012345678901234567890)
Out[47]:
BigInt

parse(BigInt, "...")从字符串输入长整数,如

In [48]:
parse(BigInt, "12345678901234567890")
Out[48]:
12345678901234567890

浮点数

现代的浮点数一般使用Float64格式, 并在内存中遵照IEEE754标准表示。 这样, 用Float64表示的浮点数可以高效地运算。

为了获得更高精度的浮点数, 可以用BigFloat类型, 这不像BigInt那样可以有任意精度, 而是有较高精度, 精度与舍入方式可以用函数setprecision()setrounding()全局性地设置。

数组处理

科学计算, 包括统计建模、数据处理, 广泛使用数组、向量、矩阵运算, 计算中大量的时间都花费在对数组的处理上面。 需要构造高效的数组处理算法。

数组存储

数组通过指定元素类型为抽象类型可以保存不同类型的数组元素, 但是这样的数组就无法利用类型信息进行优化, 所以数组应尽可能保存同一种实在类型(concrete type)的元素。

数组类型为Array{T, N}, 其中T为元素类型, 一般是实在类型, N是维数。 比如Array{Float64, 2}表示浮点数矩阵。 这种数组可以高效地存储与处理。 实在类型的数组一般可以在连续的内存中存储, 使得其遍历或者随机访问都很便捷。 而其它的动态语言如Python一般需要使用指针指向每个数组元素, 除特定的定制数组类型以外。 连续存储的数组比用指针访问的数组效率要高得多。

Julia的多维数组当元素是基础类型时按照列优先原则保存, 比如矩阵的元素也在连续的内存中保存, 先保存第一列,再保存第二列,以此类推。 多维数组类似。Fortran、Matlab也遵循列优先原则, 而C语言在保存二维数组时遵循行优先原则。

在遍历多维数组时按照其存储次序遍历效率更高; 遍历一个矩阵时, 按照列次序遍历效率更高,如:

In [49]:
A = [1 2 3; 4 5 6]
for j=1:size(A,2)
    for i=1:size(A,1)
        println("(", i, ", ", j, ") -> ", A[i,j])
    end
end
(1, 1) -> 1
(2, 1) -> 4
(1, 2) -> 2
(2, 2) -> 5
(1, 3) -> 3
(2, 3) -> 6

In [50]:
for j=1:size(A,2), i=1:size(A,1)
    println("(", i, ", ", j, ") -> ", A[i,j])
end
(1, 1) -> 1
(2, 1) -> 4
(1, 2) -> 2
(2, 2) -> 5
(1, 3) -> 3
(2, 3) -> 6

In [51]:
for i in eachindex(A)
    println(A[i])
end
1
4
2
5
3
6

下面的例子演示了按列次序与按行次序遍历矩阵的效率差别。 按列遍历的例子:

In [52]:
function f_ar01(A)
    s = zero(eltype(A))
    for j=1:size(A,2), i=1:size(A,1)
        s += A[i,j]^2
        A[i,j] = s
    end
    return s
end
f_ar01(rand(2,2))
A = rand(1_000,1_000)
@benchmark f_ar01(A)
Out[52]:
BenchmarkTools.Trial: 
  memory estimate:  16 bytes
  allocs estimate:  1
  --------------
  minimum time:     1.556 ms (0.00% GC)
  median time:      1.630 ms (0.00% GC)
  mean time:        1.652 ms (0.00% GC)
  maximum time:     9.588 ms (0.00% GC)
  --------------
  samples:          3017
  evals/sample:     1

按行遍历的例子:

In [53]:
function f_ar01b(A)
    s = zero(eltype(A))
    for i=1:size(A,1), j=1:size(A,2)
        s += A[i,j]^2
        A[i,j] = s
    end
    return s
end
f_ar01b(rand(2,2))
@benchmark f_ar01b(A)
Out[53]:
BenchmarkTools.Trial: 
  memory estimate:  16 bytes
  allocs estimate:  1
  --------------
  minimum time:     2.162 ms (0.00% GC)
  median time:      2.351 ms (0.00% GC)
  mean time:        2.414 ms (0.00% GC)
  maximum time:     8.803 ms (0.00% GC)
  --------------
  samples:          2063
  evals/sample:     1

在上面的例子中,按列遍历比按行遍历有明显的性能提升。

减少动态内存分配

对于在循环内反复实行的代码, 应避免重复地分配内存, 而是预先分配内存。 对于以数组为输入和输出的自定义函数, 可以将函数的输出也作为函数自变量。

如:

In [54]:
function f_ar02(x, y, n)
    local z
    for i in 1:n
        z = x + y
    end
    return z
end
f_ar02(rand(10), rand(10), 100)
@benchmark f_ar02(rand(10), rand(10), 10_000)
Out[54]:
BenchmarkTools.Trial: 
  memory estimate:  1.53 MiB
  allocs estimate:  10002
  --------------
  minimum time:     558.399 μs (0.00% GC)
  median time:      588.701 μs (0.00% GC)
  mean time:        629.451 μs (2.93% GC)
  maximum time:     2.763 ms (0.00% GC)
  --------------
  samples:          7914
  evals/sample:     1

进行了许多次内存分配。预先分配z的内存:

In [55]:
function f_ar02b(x, y, n)
    z = similar(x)
    for i in 1:n
        z .= x .+ y
    end
    return z
end
f_ar02b(rand(10), rand(10), 100)
@benchmark f_ar02b(rand(10), rand(10), 10_000)
Out[55]:
BenchmarkTools.Trial: 
  memory estimate:  480 bytes
  allocs estimate:  3
  --------------
  minimum time:     104.299 μs (0.00% GC)
  median time:      108.800 μs (0.00% GC)
  mean time:        110.891 μs (0.00% GC)
  maximum time:     253.101 μs (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1

后一个版本基本没有动态分配内存, 效率提升几倍。

用数组视图代替数组切片

A为矩阵, A[:,j]表示A的第j列组成的一维数组(向量), 这种子集访问称为数组切片(slice)。 对向量xx[1:end-1]这样的子集也会制作副本(切片)。

切片的问题是, A[:,j]会制作第j列的一个副本, 如果这个副本被反复多次访问, 而且未复制之前的内容在内存中没有连续存储, 复制的做法可以利用复制后连续存储快速访问的优点使得程序效率较高。 但是, 如果副本仅一次性使用, 或者数据子集本来已经是连续存储的, 复制操作会造成严重的性能损失。

例如,如下的程序对矩阵求列和:

In [56]:
function f_ar03(A)
    nc = size(A,2)
    map(j -> sum(A[:,j]), 1:nc)
end
f_ar03(rand(10,10))
@benchmark f_ar03(rand(1000,1000))
Out[56]:
BenchmarkTools.Trial: 
  memory estimate:  15.39 MiB
  allocs estimate:  1005
  --------------
  minimum time:     5.749 ms (0.00% GC)
  median time:      7.468 ms (0.00% GC)
  mean time:        8.088 ms (15.78% GC)
  maximum time:     18.919 ms (34.97% GC)
  --------------
  samples:          618
  evals/sample:     1

上面的程序中用了许多次内存分配。

在访问矩阵的一列时最好不制作副本而直接访问矩阵存储的该列。 在用到矩阵切片的代码前面冠以@views宏就可以达到此目的。如

In [57]:
function f_ar03b(A)
    nc = size(A,2)
    @views map(j -> sum(A[:,j]), 1:nc)
end
f_ar03b(rand(10,10))
@benchmark f_ar03b(rand(1000,1000))
Out[57]:
BenchmarkTools.Trial: 
  memory estimate:  7.68 MiB
  allocs estimate:  1005
  --------------
  minimum time:     4.634 ms (0.00% GC)
  median time:      5.073 ms (0.00% GC)
  mean time:        6.123 ms (14.53% GC)
  maximum time:     14.896 ms (49.12% GC)
  --------------
  samples:          816
  evals/sample:     1

第二个版本的运行效率明显提升, 内存分配减半。

SIMD并行优化

在实行类似z[i] = x[i] + y[i]这样的向量对应元素运算时, 现代CPU提供了并行实行的SIMD指令, Julia可以自动分析代码对某些代码实行SIMD优化。 可以用@simd宏要求进行SIMD优化。 这样优化的循环必须是循环各步独立的, 循环体没有分支操作,循环次数固定, 下标是连续变化的。 如:

In [58]:
function f_ar04!(z, x, y)
    n = length(x)
    for i=1:n
        z[i] = x[i] + y[i]
    end
    return nothing
end
z = zeros(10)
f_ar04!(z, rand(10), rand(10))
z = zeros(1_000_000)
@benchmark f_ar04!(z, rand(1_000_000), rand(1_000_000))
Out[58]:
BenchmarkTools.Trial: 
  memory estimate:  15.26 MiB
  allocs estimate:  4
  --------------
  minimum time:     9.780 ms (0.00% GC)
  median time:      10.554 ms (0.00% GC)
  mean time:        12.369 ms (13.84% GC)
  maximum time:     22.918 ms (36.66% GC)
  --------------
  samples:          404
  evals/sample:     1

SIMD优化版本:

In [59]:
function f_ar04b!(z, x, y)
    n = length(x)
    @inbounds @simd for i=1:n
        z[i] = x[i] + y[i]
    end
    return nothing
end
z = zeros(10)
f_ar04b!(z, rand(10), rand(10))
z = zeros(1_000_000)
@benchmark f_ar04b!(z, rand(1_000_000), rand(1_000_000))
Out[59]:
BenchmarkTools.Trial: 
  memory estimate:  15.26 MiB
  allocs estimate:  4
  --------------
  minimum time:     9.567 ms (0.00% GC)
  median time:      10.158 ms (0.00% GC)
  mean time:        12.029 ms (13.94% GC)
  maximum time:     25.226 ms (34.29% GC)
  --------------
  samples:          416
  evals/sample:     1

性能没有变化,可能是编译器本身已经包含了足够的优化。

数组遍历

在编写以数组为输入的函数时, 实际输入的数组可能是稠密数组、稀疏数组、数组的视图等。 这样,用单个下标遍历有可能降低性能。 数组一般提供eachindex()函数可以用来高效地线性遍历数组元素。

并行计算

科学计算和统计建模、数据分析中涉及大量的密集计算, 利用CPU多核心、多CPU、多台电脑组成集群并行计算是解决大规模计算问题的主要方法。

Julia语言是最近几年才发明的语言, 它吸取了现代最新的计算机科学技术成果, 语言中直接提供了对并行计算的支撑。 因为并行计算与硬件密切相关, 所以这些功能在不同运行环境中有些能发挥作用而有些受限。 目前在非用户干预的情况下Julia还是单线程计算, 例外是输入输出是非同步的, 一些内建的计算库如OpenBlas能够利用多线程。 为了实现并行计算需要用户自己调用并行计算功能。

Julia并行计算的说明见: Julia统计计算编程

亚式期权计算例子

下面的例子是用亚式Julia随机模拟方法对亚式期权定价的例子, 来自Mastering Julia。 用了不同的编程方法比较效率, Julia用了8.7秒, 未向量化的R程序用了93秒, 向量化的R程序用了12.6秒, 用Rcpp连接C++程序的版本用了7.8秒。 Julia与C++程序效率相近。

mean(x) = sum(x)/length(x)
function run_asian(N = 100000, PutCall='C';)
  # 来自Mastering Juia P.16
  # 股票期权分为看涨期权(call option)和看跌期权(put option)。
  # 交易双方称为granter(期权卖家)和beneficiary(期权买家)。
  # 以看涨期权为例,期权买家付出一定的期权价格从期权买家手里购买一个权力,
  # 在指定的时间以后可以用预先指定价格从期权卖家手里购买股票,
  # 当然对买家不利时买家可以不行使此权利。
  # 亚式期权的收益依赖于整个价格路径的算数平均,所以需要模拟许多条路径
  println("Setting option parameters")
  S0 = 100  # spot price,当前价格
  K = 100   # strike price,实行价格
  r = 0.05  # risk free rate,无风险利率
  q = 0.0   # dividend yield,红利
  v = 0.2   # Volatility,波动率
  tma = 0.25 # time to maturity,期权的交割时间

  Averaging = 'A' # 'A' 用算数平均,'G'用几何平均
  OptType = (PutCall == 'C' ? "CALL" : "PUT") # 看涨期权还是看跌期权
  println("Option type is $OptType")
  # 模拟设置
  println("Setting simulation parameters")
  T = 100   # number of time steps,模拟路径时的时间片段个数
  dt = tma / T  # time increment per step,模拟路径时每隔多长时间更新状态

  # Initialize the terminal stock price matrices
  # for the Euler and Milstein discretization schemes.
  # 初始化终端市场股票价格矩阵,采用Euler和Milstein离散化方案。
  # 矩阵S的每行为一条轨道,每列为一个时间点
  S = zeros(Float64, N, T)
  for n=1:N
    S[n,1] = S0
  end

  # Simulate the stock price under the Euler and Milstein schemes.
  # Take average of terminal stock price.
  # 下面对每条轨道(变量n),沿时间t向前模拟。最后计算终端股价的算数平均值。
  println("Looping $N times")
  A = zeros(Float64, N); # A用来存储每条轨道沿时间的平均值。
  for n = 1:N
    for t = 2:T
      dW = (randn(1)[1])*sqrt(dt) # N(0, \sqrt{dt})随机变量
      z0 = (r - q - 0.5*v*v)*S[n,t-1]*dt
      z1 = v*S[n, t-1]*dW
      z2 = 0.5*v*v*S[n,t-1]*dW*dW
      S[n,t] = S[n, t-1] + z0 + z1 + z2
    end
    if cmp(Averaging, 'A') == 0
      A[n] = mean(S[n, :])
    elseif cmp(Averaging, 'G') == 0
      A[n] = exp(mean(log(S[n,:])))
    end
  end

  # define the payoff(计算每条轨道的收益)
  P = zeros(Float64, N)
  if cmp(PutCall, 'C') == 0
    for n = 1:N
      P[n] = max(A[n] - K, 0);
    end
  elseif cmp(PutCall, 'P') == 0
    for n = 1:N
      P[n] = max(K - A[n], 0);
    end
  end

  # Calculate the price of the Asian option(计算亚式期权的价格)
  AsianPrice = exp(-r*tma)*mean(P)
  println("Price: ",  AsianPrice)
end
@time run_asian(1000000, 'C')
## Setting simulation parameters
## Looping 1000000 times
## Price: 2.5879421759595704
##   8.717699 seconds (100.18 M allocations: 10.455 GiB, 19.93% gc time)

Julia版本用了8.7秒。速度很好。

将Julia版本改成R版本,没有过多的向量化处理,用了93秒, 与Julia版本运行速度相差10倍以上:

run_asian <- function(N = 100000, PutCall='C'){
  # 来自Mastering Juia P.16,改变自Julia版本
  # 股票期权分为看涨期权(call option)和看跌期权(put option)。
  # 交易双方称为granter(期权卖家)和beneficiary(期权买家)。
  # 以看涨期权为例,期权买家付出一定的期权价格从期权买家手里购买一个权力,
  # 在指定的时间以后可以用预先指定价格从期权卖家手里购买股票,
  # 当然对买家不利时买家可以不行使此权利。
  # 亚式期权的收益依赖于整个价格路径的算数平均,所以需要模拟许多条路径
  cat("Setting option parameters\n")
  S0 = 100  # spot price,当前价格
  K = 100   # strike price,实行价格
  r = 0.05  # risk free rate,无风险利率
  q = 0.0   # dividend yield,红利
  v = 0.2   # Volatility,波动率
  tma = 0.25 # time to maturity,期权的交割时间

  Averaging = 'A' # 'A' 用算数平均,'G'用几何平均
  OptType = ifelse(PutCall == 'C', "CALL", "PUT") # 看涨期权还是看跌期权
  cat("Option type is",  OptType, "\n")

  # 模拟设置
  cat("Setting simulation parameters\n")
  T = 100   # number of time steps,模拟路径时的时间片段个数
  dt = tma / T  # time increment per step,模拟路径时每隔多长时间更新状态

  # Initialize the terminal stock price matrices
  # for the Euler and Milstein discretization schemes.
  # 初始化终端市场股票价格矩阵,采用Euler和Milstein离散化方案。
  # 矩阵S的每行为一条轨道,每列为一个时间点
  S = matrix(0.0, nrow=N, ncol=T)
  S[,1] = S0 # 第一列是初始时刻的价格

  # Simulate the stock price under the Euler and Milstein schemes.
  # Take average of terminal stock price.
  # 下面对每条轨道(变量n),沿时间t向前模拟。最后计算终端股价的算数平均值。
  cat("Looping",  N, "times\n")
  A = numeric(N)  # A用来存储每条轨道沿时间的平均值。
  for(n in 1:N){
    dW_arr <- rnorm(T-1, 0, sd=sqrt(dt))
    for(t in 2:T){
      dW = dW_arr[t-1] # N(0, \sqrt{dt})随机变量
      z0 = (r - q - 0.5*v*v)*S[n,t-1]*dt
      z1 = v*S[n, t-1]*dW
      z2 = 0.5*v*v*S[n,t-1]*dW*dW
      S[n,t] = S[n, t-1] + z0 + z1 + z2
    } # for t
    if(Averaging == 'A'){
      A[n] = mean(S[n,])
    } else if(Averaging == 'G'){
      A[n] = exp(mean(log(S[n,])))
    } else {
      stop("Averaging selection known!")
    } # if
  } # for n

  # define the payoff(计算每条轨道的收益)
  P = numeric(N)
  if(PutCall == 'C'){
    P = pmax(A - K, 0)
  } else if(PutCall == 'C'){
    P = pmax(K - A, 0)
  }

  # Calculate the price of the Asian option(计算亚式期权的价格)
  AsianPrice = exp(-r*tma)*mean(P);
  return(AsianPrice);
}
print(system.time(run_asian(1E6, 'C'))[3])

将一百万次重复模拟向量化, 关于时间t作循环,结果运行时间为12.6秒, 仅比Julia版本慢了不到一倍。程序如下:

run_asian_vec <- function(N = 100000, PutCall='C'){
  cat("Setting option parameters\n")
  S0 = 100  # spot price,当前价格
  K = 100   # strike price,实行价格
  r = 0.05  # risk free rate,无风险利率
  q = 0.0   # dividend yield,红利
  v = 0.2   # Volatility,波动率
  tma = 0.25 # time to maturity,期权的交割时间

  Averaging = 'A' # 'A' 用算数平均,'G'用几何平均
  OptType = ifelse(PutCall == 'C', "CALL", "PUT") # 看涨期权还是看跌期权
  cat("Option type is",  OptType, "\n")

  # 模拟设置
  cat("Setting simulation parameters\n")
  T = 100   # number of time steps,模拟路径时的时间片段个数
  dt = tma / T  # time increment per step,模拟路径时每隔多长时间更新状态

  # Initialize the terminal stock price matrices
  # for the Euler and Milstein discretization schemes.
  # 初始化终端市场股票价格矩阵,采用Euler和Milstein离散化方案。
  # 矩阵S的每行为一条轨道,每列为一个时间点
  S = matrix(0.0, nrow=N, ncol=T)
  S[,1] = S0 # 第一列是初始时刻的价格

  # Simulate the stock price under the Euler and Milstein schemes.
  # Take average of terminal stock price.
  # 下面对每条轨道(变量n),沿时间t向前模拟。最后计算终端股价的算数平均值。
  cat("Looping",  N, "times\n")
  dW_arr <- rnorm(T-1, 0, sd=sqrt(dt))
  for(t in 2:T){
    dW = rnorm(N, 0, sd=sqrt(dt))
    z0 = (r - q - 0.5*v*v)*S[,t-1]*dt
    z1 = v*S[, t-1]*dW
    z2 = 0.5*v*v*S[,t-1]*dW*dW
    S[,t] = S[, t-1] + z0 + z1 + z2
  } # for t

  if(Averaging == 'A'){
    A = rowMeans(S)
  } else if(Averaging == 'G'){
    A = exp(rowMeans(log(S)))
  } else {
    stop("Averaging selection known!")
  } # if

  # define the payoff(计算每条轨道的收益)
  if(PutCall == 'C'){
    P = pmax(A - K, 0)
  } else if(PutCall == 'C'){
    P = pmax(K - A, 0)
  }

  # Calculate the price of the Asian option(计算亚式期权的价格)
  AsianPrice = exp(-r*tma)*mean(P)
  return(AsianPrice)
}
system.time(print(run_asian_vec(1E6, 'C')))[3]

在R中用Rcpp将改写为C++的程序连接到R中计算,用了7.8秒,与Julia版本运行速度持平。 设C++源文件名为asian_option.cpp:

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
double run_asian_Cpp(int N = 100000) {
  double S0 = 100;  // spot price,当前价格
  double K = 100;   // strike price,实行价格
  double r = 0.05;  // risk free rate,无风险利率
  double q = 0.0;   // dividend yield,红利
  double v = 0.2;   // Volatility,波动率
  double tma = 0.25; // time to maturity,期权的交割时间
  char PutCall = 'C';
  RNGScope scope; // 设置随机数种子

  char Averaging = 'A'; // 'A' 用算数平均,'G'用几何平均

  // 模拟设置
  int T = 100;   // number of time steps,模拟路径时的时间片段个数
  double dt = tma / T;  // time increment per step,模拟路径时每隔多长时间更新状态

  // Initialize the terminal stock price matrices
  // for the Euler and Milstein discretization schemes.
  // 初始化终端市场股票价格矩阵,采用Euler和Milstein离散化方案。
  // 矩阵S的每行为一条轨道,每列为一个时间点
  NumericMatrix S(N, T);
  for(int i=0; i<N; i++)
    S(i,0) = S0; // 第一列是初始时刻的价格

  // Simulate the stock price under the Euler and Milstein schemes.
  // Take average of terminal stock price.
  // 下面对每条轨道(变量n),沿时间t向前模拟。最后计算终端股价的算数平均值。
  NumericVector A(N);  // A用来存储每条轨道沿时间的平均值。
  NumericVector dW_arr(T-1);
  double dW, z0, z1, z2;
  double An;
  for(int n=0; n<N; n++){
    dW_arr = rnorm(T-1, 0.0, sqrt(dt));
    for(int t=1; t < T; t++){
      dW = dW_arr[t-1]; // N(0, \sqrt{dt})随机变量
      z0 = (r - q - 0.5*v*v)*S(n,t-1)*dt;
      z1 = v*S(n,t-1)*dW;
      z2 = 0.5*v*v*S(n,t-1)*dW*dW;
      S(n,t) = S(n,t-1) + z0 + z1 + z2;
    } // for t
    if(Averaging == 'A'){
      An = 0.0;
      for(int i=0; i<T; i++){
        An += S(n,i);
      }
      A[n] = An / T;
    } else if(Averaging == 'G'){
      An = 0.0;
      for(int i=0; i<T; i++){
        An += log(S(n,i));
      }
      An /= T;
      A[n] = exp(An);
    } else {
      Rf_error("Error Averaging!");
    } // if
  } // for n

  // define the payoff(计算每条轨道的收益)
  NumericVector P(N);
  if(PutCall == 'C'){
    for(int i=0; i<N; i++){
      double Pi = A[i] - K;
      if(Pi < 0.0) Pi = 0.0;
      P[i] = Pi;
    }
  } else if(PutCall == 'C'){
    for(int i=0; i<N; i++){
      double Pi = K - A[i];
      if(Pi < 0.0) Pi = 0.0;
      P[i] = Pi;
    }
  }

  // Calculate the price of the Asian option(计算亚式期权的价格)
  An = 0.0;
  for(int n=0; n<N; n++){
    An += P[n];
  }
  An /= N;
  double AsianPrice = exp(-r*tma)*An;
  return(AsianPrice);
}

R接口:

sourceCpp(file="asian_option.cpp")
system.time(print(run_asian_Cpp(1E4)))[3]
## [1] 2.587442
## elapsed 
##    7.8
XML 地图 | Sitemap 地图