跑的时间测试有点不稳定,请问原因

请问大家有遇到过Julia跑出来的结果不稳定,就是方差比较大的问题吗?比如我重复下面的100次试验,记录每次所用的时间,结果显示总有那么几次的时间是比平均值高出非常多,请问有大佬知道是什么原因吗?怎么才能降低方差呢?

n=10^6;
rlevel = [[0;1;2;3;4;5;6;7;8;9] .//10; [99//100, 999//1000]];
klevel = [[1; 10; 100; 500] .// 10^4; [1;2;3;4;5;6;7;8;9] .// 10]

res = zeros(100,2);
Random.seed!(n)
for i in 1:100
    x0 = rand(Float64, n);
    x0sort = sort(x0, rev=true);
    k = Int64(ceil(n * klevel[1])) # 10
    tks = sum(x0sort[1:k]);
    r = tks * float(rlevel[1]); # 1

    res[i,1] = @elapsed begin
        tol::Float64 = 1e-8 # stopping criteria
        uL::Float64 = -Inf
        diffL::Float64 = vlengthL::Float64 = uR::Float64 = +Inf
        l_fold = unew::Float64 = r/k;
        uold = max_a::Float64 = maximum(x0);
        Ite_init::Integer = 0;
        sum_a_r1_old = 0.0;
        m_old = 0;
        l_fnew = 0.0;
        flag = -2;
      while true
        r1 = x0 .> unew
        m = sum(r1)
        sum_a_r1 = sum(x0[r1])
        l_fnew = (r - sum_a_r1 + m*unew)/k
    
        if unew - l_fnew < tol
            flag = 0;
            break
        end
  
        if m == k
          flag = -1;
          break
        end
        
        if m > k
          break
        end
        
        Ite_init += 1;
        uold = unew;
        l_fold = l_fnew
        m_old = m;
        sum_a_r1_old = sum_a_r1
        
        unew = (uold*m - k*l_fold)/(m-k);

      end
    end

    res[i, 2] = @elapsed begin
        r3 = x0.>l_fold
        compare_result = (r-k*l_fold) - (sum(x0[r3]) - sum(r3)*l_fold) + k*(uold - l_fold)
        flag = -1;
    end
end
mean(res, dims=1)
std(res, dims=1)

plot出来的结果如下:

感觉是自动gc, 这个是正常的, 没有特殊要求没必要降低方差

我也发现了是gc的问题。但是方差也是衡量我算法表现的一个因素,所以我就想把方差降下去。我发现了有可能是循环里面

r1 = x0 .> unew
m = sum(r1)
sum_a_r1 = sum(x0[r1])

这个的问题,产生了太多额外分配,后面改成一个for循环遍历元素好像就把方差降下去了。但是这个问题好像在matlab和python里面是没有的。

几个常见的点:

  1. 写算法测试最好封装在函数里面。在global scope里面大量的优化是不会进行的。
  2. 检查函数类型稳定(@code_warntype)。
  3. @time 或者其他工具查看allocation以及GC timing。尽量做到内部循环里面不要有新的allocation。

谢谢你的建议。

  1. 其实我已经封装在函数里面了,但是在做实验中发现了不稳定的问题,就把算法的一部分拿出来去检测,就是我放的代码。
  2. 已经用了这个@code_warntype,但是并没有显示什么问题。
  3. 是的,现在我正在减少新的allocation。
1 个赞

如果你怀疑的GC是导致时间测试变化很大的原因,在@time 里面有一项time spent in GC估计会对应显示很大的时间比例用在了GC上。可以看看是不是这样。

part1函数就是我提供代码里面的while循环,封装成的函数。

如果while循环里面使用原来的这个代码:

r1 = x0 .> unew
m = sum(r1)
sum_a_r1 = sum(x0[r1])

那么得到的结果是

julia> @time part1(x0, r, k)
  0.054270 seconds (114.02 k allocations: 8.540 MiB, 17.76% gc time, 99.40% compilation time)

如果把上述的代码用一个for循环遍历,

m = 0
sum_a_r1 = 0.0
@inbounds @simd for xi in x0
    if xi > unew
    m += 1
    sum_a_r1 += xi
    end
end

那么得到的结果是:

julia> @time part1(x0, r, k)
  0.011067 seconds (10.11 k allocations: 708.766 KiB, 98.65% compilation time)

感觉就是在循环里面产生了额外的分配导致的这个时间区别

你的第一种写法在Python和MATLAB里面做vectorization非常常见,但在Julia里面的确应该写成更接近C和Fortran的方式。但另外,Julia里面的sum其实有更高级的用法:

sum_a_r1 = sum(x -> x > unew, x0)

谢谢!我对比了一下sum这个写法和我用for循环求sum的写法,时间上面是差不多的。