显式定义变量与访问数组的区别

Question

显式定义变量与访问数组的区别

3

我正在实现具有自适应步长（RK45）的龙格-库塔-费尔伯格方法。我在笔记本中定义并调用了Butcher表格。

module FehlbergTableau
    using StaticArrays
    export A, B, CH, CT

    A = @SVector [ 0 , 2/9 , 1/3 , 3/4 , 1 , 5/6 ]

    B = @SMatrix [ 0          0         0        0        0
                   2/9        0         0        0        0
                   1/12       1/4       0        0        0
                   69/128    -243/128   135/64   0        0
                  -17/12      27/4     -27/5     16/15    0
                   65/432    -5/16      13/16    4/27     5/144 ]

    CH = @SVector [ 47/450 , 0 , 12/25 , 32/225 , 1/30 , 6/25 ]

    CT = @SVector [ -1/150 , 0 , 3/100 , -16/75 , -1/20 , 6/25 ]
end

using .FehlbergTableau

如果我按照 RK45 算法的直接编码方式实现它：

function infinitesimal_flow(A::SVector{6,Float64}, B::SMatrix{6,5,Float64}, CH::SVector{6,Float64}, CT::SVector{6,Float64}, t0::Float64, Δt::Float64, J∇H::Function, x0::SVector{N,Float64}) where N
    
    k1 = Δt * J∇H( t0 + Δt*A[1], x0 ) 
    k2 = Δt * J∇H( t0 + Δt*A[2], x0 + B[2,1]*k1 ) 
    k3 = Δt * J∇H( t0 + Δt*A[3], x0 + B[3,1]*k1 + B[3,2]*k2 )
    k4 = Δt * J∇H( t0 + Δt*A[4], x0 + B[4,1]*k1 + B[4,2]*k2 + B[4,3]*k3 )
    k5 = Δt * J∇H( t0 + Δt*A[5], x0 + B[5,1]*k1 + B[5,2]*k2 + B[5,3]*k3 + B[5,4]*k4 )
    k6 = Δt * J∇H( t0 + Δt*A[6], x0 + B[6,1]*k1 + B[6,2]*k2 + B[6,3]*k3 + B[6,4]*k4 + B[6,5]*k5 )
    
    TE = CT[1]*k1 + CT[2]*k2 + CT[3]*k3 + CT[4]*k4 + CT[5]*k5 + CT[6]*k6  
    xt = x0 + CH[1]*k1 + CH[2]*k2 + CH[3]*k3 + CH[4]*k4 + CH[5]*k5 + CH[6]*k6 
    
    norm(TE), xt
end

并将其与更紧凑的实现进行比较

function infinitesimal_flow_2(A::SVector{6,Float64}, B::SMatrix{6,5,Float64}, CH::SVector{6,Float64}, CT::SVector{6,Float64}, t0::Float64,Δt::Float64,J∇H::Function, x0::SVector{N,Float64}) where N
    
    k = MMatrix{N,6}(0.0I)
    TE = zero(x0); xt = x0
    
    for i=1:6
        # EDIT: this is wrong! there should be a new variable here, as pointed 
        # out by Lutz Lehmann: xs = x0
        for j=1:i-1
            # xs += B[i,j] * k[:,j]
            x0 += B[i,j] * k[:,j] #wrong
        end
        k[:,i] = Δt *  J∇H(t0 + Δt*A[i], x0) 

        TE += CT[i]*k[:,i]
        xt += CH[i]*k[:,i]B[i,j] * k[:,j]
    end
    norm(TE), xt
end

那么第一个显式定义变量的函数会更快：

J∇H(t::Float64, X::SVector{N,Float64}) where N = @SVector [ -X[2]^2, X[1] ]
x0 = SVector{2}([0.0, 1.0])
infinitesimal_flow(A, B, CH, CT, 0.0, 1e-2, J∇H, x0)
infinitesimal_flow_2(A, B, CH, CT, 0.0, 1e-2, J∇H, x0)

@btime infinitesimal_flow($A, $B, $CH, $CT, 0.0, 1e-2, $J∇H, $x0)
>> 19.387 ns (0 allocations: 0 bytes)
@btime infinitesimal_flow_2($A, $B, $CH, $CT, 0.0, 1e-2, $J∇H, $x0)
>> 50.985 ns (0 allocations: 0 bytes)

我无法找到类型不稳定或任何可以证明滞后的理由，对于更复杂的表格，必须使用循环形式的算法。我做错了什么？

附言：infinitesimal_flow_2中的瓶颈是k[:,i] = Δt * J∇H(t0 + Δt*A[i], x0)这一行。

- QuantumBrick

一个2.63倍的减速相对来说还算温和。你的@btime结果表明编译器消除了MMatrix及其切片的分配，而且J∇H在两个函数中以相同的参数类型运行6次，所以我不认为这是原因。我能想到的唯一剩下的事情就是for循环增加了一点开销，并使编译器难以优化（就像它可以在展开版本中识别xt = x0 + TE一样），但这太小了，我不能确定它是否会使运行时间减少2.63倍。如果只是阅读和比较@code_native那该多好啊！ - BatWannaBe

你所说的自由度具体指什么？是N吗？6变成12了吗？如果它增加了你的"MMatrix"大小，那可能就是为什么它变慢1000倍的原因；分配消除总是有点神秘的，而且编译器可能会停止对更大尺寸进行优化。也许你应该发表12个自由度版本及其"@btime"结果？"1000x"的减速度会更引人注目。 - BatWannaBe

我不确定你为什么说它是不可变的。MMatrix 是一种可变类型，且在堆上分配，文档也是这样说的。你还在这一行中对其进行了修改：k[:,i] = Δt * J∇H(t0 + Δt*A[i], x0)。我的观点只是堆分配是减慢速度的一个重要来源，并且使可变数组 k 和所有复制的切片 k[:,i] 通常会分配一堆内存。编译器在当前帖子的示例中做了相当多的工作，将其减少到“0次分配”（意味着在堆上）。 - BatWannaBe

@BatWannaBe 你是对的，矩阵确实会变大。只不过这不是因为6（这个索引是固定的），而是因为N。我现在同意你的观点。 - QuantumBrick

1

另外，从展开版本来看，似乎您不需要矩阵k的零值。可以这样定义：k = MMatrix{N,6, Float64, 6*N}(undef)。虽然在您的示例中看起来并不多，但对于更大的矩阵，收益更为显著。 - Thomas Jalabert

显示剩余4条评论

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Lutz Lehmann · Accepted Answer

每个RK方法的阶段都直接从RK步骤的基础点计算其评估点。在第一种方法中，这是显式的。在第二种方法中，您需要在每个阶段重新设置点计算，例如在...

    for i=1:6
        xs = x0
        for j=1:i-1
            xs += B[i,j] * k[:,j]
        end
        k[:,i] = Δt *  J∇H(t0 + Δt*A[i], xs) 
        ...

最小的步骤计算错误都可能会严重影响步长控制器，导致步长逐渐减小并使得计算量急剧增加。一个例子是RKF45中的4101错误。