Juliaで多次元配列の導関数を計算する効率的な方法を探しています。正確には、私は同等のものが欲しいです numpy.gradient
ジュリアで。ただし、ジュリア関数 diff
:
- 2次元配列でのみ機能します
- 微分次元に沿って配列のサイズを1つ減らします
の定義を拡張するのは簡単です diff
ジュリアの3次元配列で動作できるように、例えば。と
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
これは、たとえば
a = [i*j*k for i in 1:10, j in 1:10, k in 1:20]
ただし、任意の次元への拡張は不可能であり、境界は考慮されないため、グラデーションは元の配列と同じ次元を持つことができます。
私はのアナログを実装するためのいくつかのアイデアがありますジュリアでのnumpyの勾配ですが、非常に遅くて醜いのではないかと心配しています。したがって、私の質問:ジュリアでこれを行うための標準的な方法はありますか?
ありがとう。
回答:
回答№1は4あまりに慣れていない diff
、しかし、その実行内容について私が理解していることから、パラメトリックタイプやスプラッティングなどのJulia機能を使用するn次元の実装を作成しました。
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
いくつかの健全性チェックを伴う:
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.
コード生成を使用して有限次元まで特殊なメソッドを作成するなど、これを行うためのより魔法の方法があることに注意してください。ただし、十分なパフォーマンスを得るには、おそらくそれは必要ないと思います。
回答№2の場合は1
それを行うさらに簡単な方法:
mydiff(A::AbstractArray,dim) = mapslices(diff, A, dim)
しかし、これが速度の点でどのように比較されるかはわかりません。
編集: 少し遅いかもしれませんが、これは関数を高階配列に拡張するためのより一般的な解決策です。
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)