Hello, 
I need to calculate the gradient of a 3 dimensional matrix. Is there today in 2023 an equivalent of numpy.gradient in Julia?(numpy.gradient — NumPy v1.24 Manual)
This question had already been asked, but the answer was not really satisfactory I think.
Thanks for your answer !
fdekerm
What are you missing in the answers given before?
Unless I missed something, I find that all the proposed solutions are quite complex and “artisanal” for something as simple as this, and which is done in one line with numpy. Especially since there was a time when there was a “gradient()” function directly in Base.
OK, if I understand you correctly you would like to have a simple gradient() utility function that works on any number of dimensions, and you do not know in which package it should be placed, correct?
Yes, exactly!
I did not find any related package on JuliaHub .
Perhaps it would be worth to create a new package with such a function? Any comments?
I added a gradient function to CoupledFields.jl (]add CoupledFields#master to try, see Example 3 in the docs). It uses a different method than numpy.gradient, but has the advantage of working for regular grids or irregular gridpoints. Let me know if it works for your case!
@Mattriks I’m confused by your syntax.
If I try using CoupledFields: ndgrid, I get the error
ERROR: UndefVarError: ndgrid not defined
I also get errors for i = range(-2, 2, 6). How do I run your example?
Can you please tell which exactly method do you use for gradient calculation in case of irregular grid points.
For now, the new function is available on master branch:
]add CoupledFields#master
I’m going to add a “Theory” page to the CoupledFields documentation, but for now there is a brief explanation here: Smooth delta functions to compute Jacobian and Hessian from samples - #10 by Mattriks
The fact that this is hard to find may be because in Julia you often don’t actually need to store the gradient as an array (or tuple-of-arrays), instead you compute it pointwise on the fly as an input to whatever other computation you’re actually doing. In other words, the gradient does not need to be an “endpoint” of something in Julia.
That said, there are a couple that haven’t been mentioned:
imgradients(A, KernelFactors.sobel) should basically give you what you want, but know that sobel does a bit of averaging over nearest-neighbors and will not be exactly the same as the two-point gradient you’re currently using.The OP mentions “the gradient of a 3d matrix” (which I assume means 3d positions, and a 1d variable at those positions), can imgradients be used?
Interpolations.jl works for this, too, although with a bit more overhead than in-place calculation:
julia> using Interpolations
julia> v = rand(3, 3, 3);
julia> itp = linear_interpolation(axes(v), v);
julia> [Interpolations.gradient(itp, idx...) for idx in knots(itp)]
3×3×3 Array{StaticArraysCore.SVector{3, Float64}, 3}:
[:, :, 1] =
[0.0217651, 0.081234, 0.245227] [0.422528, 0.50527, 0.237876] [-0.336426, 0.50527, -0.153063]
[-0.140606, 0.481997, 0.47315] [-0.205345, -0.253683, -0.562774] [-0.184544, -0.253683, -0.330334]
[-0.140606, 0.417258, 0.815596] [-0.205345, -0.232882, -0.0649707] [-0.184544, -0.232882, 0.0365438]
[:, :, 2] =
[0.249688, 0.0738826, 0.347969] [-0.378122, 0.114332, -0.0942331] [-0.513697, 0.114332, -0.253928]
[0.201841, -0.553927, 0.189845] [0.292458, -0.0212431, 0.726295] [0.182334, -0.0212431, 0.368639]
[0.201841, -0.463309, 0.0285594] [0.292458, -0.131368, 0.152346] [0.182334, -0.131368, 0.292275]
[:, :, 3] =
[0.0915639, -0.368319, 0.347969] [0.442407, -0.0453636, -0.0942331] [0.10887, -0.0453636, -0.253928]
[0.0405551, -0.0174767, 0.189845] [-0.281491, -0.3789, 0.726295] [0.105971, -0.3789, 0.368639]
[0.0405551, -0.339523, 0.0285594] [-0.281491, 0.00856203, 0.152346] [0.105971, 0.00856203, 0.292275]
Yes, it’s arbitrary-dimensional.
What about 3d positions and a 2d variable? 
Not fully sure I know what you mean, is it this?
julia> using StaticArrays, ImageFiltering
julia> A = [SVector(rand(), rand()) for i = 1:3, j = 1:3, k = 1:3];
julia> g1, g2, g3 = imgradients(A, KernelFactors.sobel);
julia> g1
3×3×3 Array{SVector{2, Float64}, 3}:
[:, :, 1] =
[0.0762768, 0.122968] [0.00718473, 0.081613] [-0.0166742, -0.0798793]
[0.127677, 0.16464] [0.0277456, 0.130215] [-0.154041, 0.105112]
[0.0514005, 0.0416719] [0.0205609, 0.0486023] [-0.137366, 0.184991]
[:, :, 2] =
[0.0474921, 0.155314] [0.0466046, 0.120711] [-0.0332256, 0.0171544]
[0.137826, 0.202603] [0.113671, 0.134971] [0.000161656, 0.0273804]
[0.0903337, 0.0472892] [0.0670666, 0.0142602] [0.0333872, 0.010226]
[:, :, 3] =
[0.00795333, -0.0351133] [0.066914, -0.0367565] [-0.0944003, -0.0250204]
[0.0625429, 0.102988] [0.112584, 0.106004] [-0.0335448, -0.0150227]
[0.0545896, 0.138101] [0.0456702, 0.14276] [0.0608555, 0.00999772]
It also works with A = [rand(SMatrix{5,5}) for i = 1:3, j = 1:3, k = 1:3];, if that’s what you mean by a 2d variable.
That looks right, I should do a comparison sometime.
Keep in mind that you may see small numerical differences due to the fact that this notion of “gradient” comes from applying (in 2d) the following stencils:
julia> using ImageFiltering
julia> g1, g2 = Kernel.sobel();
julia> g1
3×3 OffsetArray(::Matrix{Float64}, -1:1, -1:1) with eltype Float64 with indices -1:1×-1:1:
-0.125 -0.25 -0.125
0.0 0.0 0.0
0.125 0.25 0.125
julia> g2
3×3 OffsetArray(::Matrix{Float64}, -1:1, -1:1) with eltype Float64 with indices -1:1×-1:1:
-0.125 0.0 0.125
-0.25 0.0 0.25
-0.125 0.0 0.125
You should be able to replicate the numpy notion exactly with
g1, g2 = KernelFactors.sobel((true, false), 1), KernelFactors.sobel((false, true), 2)
which is just a (too complicated?) way of imfiltering with centered(SVector(-0.5, 0, 0.5)) and centered(SVector(-0.5, 0, 0.5)') (which doesn’t quite work).
It’d be great to have these standard finite difference kernels in ImageFiltering or some extension package, and/or to spin out the kernel/window functions into a separate package (KernelOperators.jl, maybe under the JuliaArrays org?). At the moment, I don’t think there’s a package for textbook vector calculus operators on gridded data, and I’ve ended up cobbling them together whenever they’re needed from Interpolations or ImageFiltering. For curl in particular, this requires wasteful allocation/computation of the diagonal part of the Jacobian.