assegnazioni non necessari utilizzando gli operatori aggiornamento Julia

Question

Si consideri la seguente funzione:

function mytest(x, b) 
    y = zeros(x[:,:,1]) 
    for i in 1:length(b) 
     y += b[i] * x[:,:,i] 
    end 
    return y 
end 

quando l'eseguo, ricevo il seguente:

x = rand(30,30,100000) 
b = rand(100000) 
@time mytest(x,b) 

elapsed time: 0.571765222 seconds (727837732 bytes allocated, 66.49% gc time) 

Perché è allocare la memoria così tanto e spendere così tanto tempo fare la raccolta dei rifiuti? Il codice dovrebbe essere di tipo stabile e mi aspetto che l'operatore += non esegua alcuna ridistribuzione. Tuttavia, sembra che si stia ridistribuendo ogni volta che aggiunge due matrici.

Devo considerare che si tratta di un bug in Julia? E ancora più importante, come posso scrivere questo codice in un modo che non si riassegna?

MODIFICA: corretto errore di battitura.

Answer 1

@ cd98 ha richiesto la mia soluzione a ciclo triplo, che risolve il problema di allocazione, ma presumo che sottolinei una versione vettoriale equivalente. Eccolo:

function mytest(x, b) 
    d1, d2, d3 = size(x) 
    y = zeros(eltype(x), d1, d2) 
    for i in 1:d3 
     for j in 1:d2 
      for k in 1:d1 
       y[k,j] += b[i] * x[k,j,i] 
      end 
     end 
    end 
    return y 
end 

x = rand(30,30,100000) 
b = rand(100000) 
@time mytest(x,b) 
@time mytest(x,b) 

E l'output:

elapsed time: 0.218220119 seconds (767172 bytes allocated) 
elapsed time: 0.197181799 seconds (7400 bytes allocated) 

Answer 2

V'è inoltre la possibilità di utilizzare Base.Cartesian: dopo using Base.Cartesian, è possibile scrivere

@nloops 3 i x begin 
    (@nref 2 y i) += b[i_3] * (@nref 3 x i) 
end 

che si espande per fondamentalmente la stessa loop come nella risposta di Jim.

Answer 3

non risolve il (originale) problema di allocazione, ma io ottenere circa un aumento di velocità 1,8 volte semplicemente avvolgendo i loop in quest'ultima soluzione @inbounds:

function mytest_inbounds(x, b) 
    d1, d2, d3 = size(x) 
    y = zeros(eltype(x), d1, d2) 
    @inbounds begin 
     for i in 1:d3 
      for j in 1:d2 
       for k in 1:d1 
        y[k,j] += b[i] * x[k,j,i] 
       end 
      end 
     end 
    end 
    return y 
end 

x = rand(30, 30, 100000) 
b = rand(100000) 
@time mytest(x, b) 
@time mytest(x, b) 
@time mytest_inbounds(x, b) 
@time mytest_inbounds(x, b) 

uscita:

elapsed time: 0.39144919 seconds (767212 bytes allocated) 
elapsed time: 0.353495867 seconds (7400 bytes allocated) 
elapsed time: 0.202614643 seconds (396972 bytes allocated) 
elapsed time: 0.193425902 seconds (7400 bytes allocated) 

Inoltre, un sacco di buona discussione sulle questioni connesse qui:

https://groups.google.com/forum/#!msg/julia-users/aYS_AvKqPCI/DyTiq4lKIAoJ