Differentiation without explicit function (np.gradient)

New to Julia
,
1

Hi everyone,
I am currently looking to evaluate some numerical data and need a method similar to numpy.gradient. All the packages, I have found seem to assume that I have some continuous function that I can evaluate at arbitrary points, however I only have a fixed array of function values that I can use (i.e. x = Array{Float64,1}, y = Array{Float64,1}).
I have found an old topic (Is there a central difference/gradient function somewhere?) that discusses this, but I found that I was not completely satisfied with the discussed solutions and hope that something has changed since then.
Most packages seem to have some advanced functions to differentiate matrices or higher dimensional arrays, but not simple 1-D arrays.
The closest answer is this quick and dirty code of mine, inspired by one of the answers from above that gives what I want for evenly spaced x -Arrays:

# returns derivative of f wrt. x such that the result has a length equal to that of x
function deriv(y::AbstractVector,x::AbstractVector)
    function centraldiff(v::AbstractVector)
        dv  = diff(v)/2 # half the derivative
        a   = [dv[1];dv] # copies first element
        a .+= [dv;dv[end]] # copies last element, add both results to compute average
        return(a)
    end
    return  centraldiff(y)./centraldiff(x)
end

but for such a simple problem I would not expect having to rely on my own clumsy implementation. It would be nice to import some function that can do this for non-evenly spaced x as well.
Have I overlooked a suitable function in one of the many differentiation libraries, or is my problem ill-formulated somehow, such that there is a much better workflow altogether? Numerically, it is not feasible to compute values of y dynamically.
I hope you can help me in this problem!

1 Like
2

One nice option is to use Interpolations.jl to create a piecewise linear interpolation and then ask for its derivatives:

julia> using Interpolations

julia> x = [1, 2, 4];

julia> y = [1, 2, 3];

julia> itp = interpolate((x,), y, Gridded(Linear()));

You can now use itp to interpolate values:

julia> itp(3.0)
2.5

and compute derivatives:

julia> Interpolations.gradient(itp, 3.0)
1-element StaticArrays.SArray{Tuple{1},Float64,1,1} with indices SOneTo(1):
 0.5

Note that this gives you a gradient vector, but you can get the scalar derivative by taking its only element:

julia> only(Interpolations.gradient(itp, 3.0))
0.5

and you can use broadcasting to get multiple derivatives at different x values:

julia> Interpolations.gradient.(Ref(itp), [1.0, 2.0, 3.0])
3-element Array{StaticArrays.SArray{Tuple{1},Float64,1,1},1}:
 [1.0]
 [1.0]
 [0.5]

or as scalars:

julia> only.(Interpolations.gradient.(Ref(itp), [1.0, 2.0, 3.0]))
3-element Array{Float64,1}:
 1.0
 1.0
 0.5
10 Likes
3

Oh, and you can even use your interpolation as if it were a continuous-valued function (because it is), so tools like ForwardDiff also just work:

julia> using ForwardDiff

julia> ForwardDiff.derivative(itp, 3.0)
0.5

although it won’t be any faster than using Interpolations.gradient.

2 Likes
4

Thank you, this is the kind of elegant solution I was looking for!

5

I think this usage of Ref() for evaluating the gradient at multiple points should be included in the documentation of Interpolations.jl. When you construct an interpolator, I think the desire to evaluate the gradient at multiple points is quite common, and yet I had no idea how to achieve it without using a for loop over the sampling points.

2 Likes
6

Note that it is not specific to that package, but a general way to escape broadcasting.

1 Like
7

Use a list comprehension like [Interpolations.gradient(itp, x) for x in [1.0, 2.0, 3.0]]

2 Likes
8

Thank you for your answer, I was looking for the same answer for a while now! The only thing is that “only” is not defined within the Interpolations package and apparently not even within the Julia language. How can this be possible? Can it be the case that it was changed after the updates?

9

According to the docs you need at least Julia 1.4, maybe that is the issue in your case?

help?> only
search: only ReadOnlyMemoryError isreadonly download Clonglong @noinline functionloc UnionAll Culonglong @functionloc Rational

  only(x)

  Returns the one and only element of collection x, and throws an ArgumentError if the collection has zero or multiple elements.

  See also: first, last.

  │ Julia 1.4
  │
  │  This method requires at least Julia 1.4.

julia> only([1])
1