Wielowymiarowy diff / gradient w Julii - tablica wielowymiarowa, julia-lang, metody numeryczne

Szukam efektywnego sposobu obliczania pochodnych macierzy wielowymiarowej w Julii. Mówiąc ściślej, chciałbym mieć równowartość numpy.gradient w Julii. Jednak funkcja Julii diff :

• działa tylko dla tablic dwuwymiarowych
• zmniejsza rozmiar tablicy o jeden wzdłuż zróżnicowanego wymiaru

Po prostu rozszerz definicję diff Julii, dzięki czemu może działać na trójwymiarowych tablicach, np. z

function diff3D(A::Array, dim::Integer)
if dim == 1
[A[i+1,j,k] - A[i,j,k] for i=1:size(A,1)-1, j=1:size(A,2), k=1:size(A,3)]
elseif dim == 2
[A[i,j+1,k] - A[i,j,k] for i=1:size(A,1), j=1:size(A,2)-1, k=1:size(A,3)]
elseif dim == 3
[A[i,j,k+1] - A[i,j,k] for i=1:size(A,1), j=1:size(A,2), k=1:size(A,3)-1]
else
throw(ArgumentError("dimension dim must be 1, 2, or 3 got $dim"))
end
end

który działałby np.

a = [i*j*k for i in 1:10, j in 1:10, k in 1:20]

Jednak rozszerzenie na dowolny wymiar nie jest możliwe, a granica nie jest brana pod uwagę, więc gradient może mieć taki sam wymiar jak pierwotna tablica.

Mam kilka pomysłów na implementację analogugradient numpy w Julii, ale obawiam się, że byłyby one bardzo powolne i brzydkie, stąd moje pytania: czy istnieje kanoniczny sposób na zrobienie tego w Julii, za którym tęskniłem? A jeśli nie, to co byłoby optymalne?

Dzięki.

Odpowiedzi:

4 dla odpowiedzi № 1

Nie jestem zbyt obeznany diff, ale z tego, co rozumiem na temat jego działania, stworzyłem n-wymiarową implementację, która wykorzystuje funkcje Julii, takie jak typy parametryczne i rozpryskiwanie:

function mydiff{T,N}(A::Array{T,N}, dim::Int)
@assert dim <= N
idxs_1 = [1:size(A,i) for i in 1:N]
idxs_2 = copy(idxs_1)
idxs_1[dim] = 1:(size(A,dim)-1)
idxs_2[dim] = 2:size(A,dim)
return A[idxs_2...] - A[idxs_1...]
end

z pewnymi sprawdzeniami czystości:

A = rand(3,3)
@assert diff(A,1) == mydiff(A,1)  # Base diff vs my impl.
@assert diff(A,2) == mydiff(A,2)  # Base diff vs my impl.

A = rand(3,3,3)
@assert diff3D(A,3) == mydiff(A,3)  # Your impl. vs my impl.

Zauważ, że są na to bardziej magiczne sposoby, takie jak generowanie kodu do tworzenia wyspecjalizowanych metod aż do skończonego wymiaru, ale myślę, że prawdopodobnie nie jest to konieczne, aby uzyskać wystarczająco dobrą wydajność.

---

1 dla odpowiedzi nr 2

Jeszcze prostszy sposób to zrobić:

mydiff(A::AbstractArray,dim) = mapslices(diff, A, dim)

Nie jestem jednak pewien, jak by to wyglądało pod względem prędkości.

Edytować: Może nieco wolniej, ale jest to bardziej ogólne rozwiązanie rozszerzające funkcje na tablice wyższego rzędu:

julia> using BenchmarkTools

julia> function mydiff{T,N}(A::Array{T,N}, dim::Int)
@assert dim <= N
idxs_1 = [1:size(A,i) for i in 1:N]
idxs_2 = copy(idxs_1)
idxs_1[dim] = 1:(size(A,dim)-1)
idxs_2[dim] = 2:size(A,dim)
return A[idxs_2...] - A[idxs_1...]
end
mydiff (generic function with 1 method)

julia> X = randn(500,500,500);

julia> @benchmark mydiff($X,3)
BenchmarkTools.Trial:
samples:          3
evals/sample:     1
time tolerance:   5.00%
memory tolerance: 1.00%
memory estimate:  2.79 gb
allocs estimate:  22
minimum time:     2.05 s (15.64% GC)
median time:      2.15 s (14.62% GC)
mean time:        2.16 s (11.05% GC)
maximum time:     2.29 s (3.61% GC)

julia> @benchmark mapslices(diff,$X,3)
BenchmarkTools.Trial:
samples:          2
evals/sample:     1
time tolerance:   5.00%
memory tolerance: 1.00%
memory estimate:  1.99 gb
allocs estimate:  3750056
minimum time:     2.52 s (7.90% GC)
median time:      2.61 s (9.17% GC)
mean time:        2.61 s (9.17% GC)
maximum time:     2.70 s (10.37% GC)

---