algorithm - 定点算法中的内存分配

Question

我需要找到函数 f 的不动点。算法很简单:

1. 给定 X，计算 f(X)
2. 如果||X-f(X)||低于一定容忍度，退出返回X，否则设置 X 等于 f(X) 并回到 1。

我想确保我不会在每次迭代时都为新对象分配内存

目前，算法如下所示:

iter1 = function(x::Vector{Float64})
    for iter in 1:max_it
        oldx = copy(x)
        g1(x)
        delta = vnormdiff(x, oldx, 2)
        if delta < tolerance
            break
        end
    end
end

这里的g1(x)是一个函数，将x设置为f(x)

但似乎这个循环在每个循环中都分配了一个新向量(见下文)。

另一种编写算法的方法如下:

iter2 = function(x::Vector{Float64})
    oldx = similar(x)
    for iter in 1:max_it
        (oldx, x) = (x, oldx)
        g2(x, oldx)
        delta = vnormdiff(oldx, x, 2)
        if delta < tolerance
            break
        end
    end
end

其中 g2(x1, x2) 是一个将 x1 设置为 f(x2) 的函数。

这是编写此类迭代问题的最有效、最自然的方法吗？

Edit1:时序显示第二个代码更快:

using NumericExtensions
max_it = 1000
tolerance = 1e-8
max_it = 100

g1 = function(x::Vector{Float64}) 
    for i in 1:length(x)
        x[i] = x[i]/2
    end
end

g2 = function(newx::Vector{Float64}, x::Vector{Float64}) 
    for i in 1:length(x)
        newx[i] = x[i]/2
    end
end

x = fill(1e7, int(1e7))
@time iter1(x)
# elapsed time: 4.688103075 seconds (4960117840 bytes allocated, 29.72% gc time)
x = fill(1e7, int(1e7))
@time iter2(x)
# elapsed time: 2.187916177 seconds (80199676 bytes allocated, 0.74% gc time)

Edit2:使用 copy!

iter3 = function(x::Vector{Float64})
    oldx = similar(x)
    for iter in 1:max_it
        copy!(oldx, x)
        g1(x)
        delta = vnormdiff(x, oldx, 2)
        if delta < tolerance
            break
        end
    end
end
x = fill(1e7, int(1e7))
@time iter3(x)
# elapsed time: 2.745350176 seconds (80008088 bytes allocated, 1.11% gc time)

Answer 1

我想替换第一段代码中的以下几行

for iter = 1:max_it
    oldx = copy( x )
    ...

通过

oldx = zeros( N )
for iter = 1:max_it
    oldx[:] = x    # or copy!( oldx, x )
    ...

会更有效率，因为没有分配数组。此外，可以通过显式编写 for 循环来提高代码的效率。例如从下面的对比可以看出这一点

function test()
    N = 1000000

    a = zeros( N )
    b = zeros( N )

    @time c = copy( a )

    @time b[:] = a

    @time copy!( b, a )

    @time \
    for i = 1:length(a)
        b[i] = a[i]
    end

    @time \
    for i in eachindex(a)
        b[i] = a[i]
    end
end

test()

在Linux(x86_64)上用Julia0.4.0得到的结果是

elapsed time: 0.003955609 seconds (7 MB allocated)
elapsed time: 0.001279142 seconds (0 bytes allocated)
elapsed time: 0.000836167 seconds (0 bytes allocated)
elapsed time: 1.19e-7 seconds (0 bytes allocated)
elapsed time: 1.28e-7 seconds (0 bytes allocated)

似乎 copy!() 比在左侧使用 [:] 更快，尽管在重复计算中差异变得微不足道(似乎有第一次 [:] 计算的一些开销)。顺便说一句，最后一个使用 eachindex() 的例子对于循环多维数组非常方便。

可以对 vnormdiff() 进行类似的比较，其中使用 norm( x - oldx ) 等比向量范数的显式循环慢，因为前者为 x - oldx 分配一个临时数组。