DiffEqDiffTools.jl is a component package in the DifferentialEquations ecosystem. It holds the common tools for taking derivatives, gradients, Jacobians, Hessians, etc. This library is for maximizing speed while giving a usable interface to end users. Included is:
If you want the fastest versions, create a cache and repeatedly call the
differencing functions at different x
values (or with different f
functions),
while if you want a quick and dirty numerical answer, directly call a differencing
function.
The general structure of the library is as follows. You can call the differencing functions directly and this will allocate a temporary cache to solve the problem with. To make this non-allocating for repeat calls, you can call the cache construction functions. Each cache construction function has two possibilities: one version where you give it prototype arrays and it generates the cache variables, and one fully non-allocating version where you give it the cache variables. This is summarized as:
In all functions, the inplace form is f!(dx,x)
while the out of place form is dx = f(x)
.
colorvec vectors are allowed to be supplied to the Jacobian routines, and these are
the directional derivatives for constructing the Jacobian. For example, an accurate
NxN tridiagonal Jacobian can be computed in just 4 f
calls by using
colorvec=repeat(1:3,N÷3)
. For information on automatically generating colorvec
vectors of sparse matrices, see SparseDiffTools.jl.
Hessian coloring support is coming soon!
DiffEqDiffTools.finite_difference_derivative(f, x::T, fdtype::Type{T1}=Val{:central},
returntype::Type{T2}=eltype(x), f_x::Union{Nothing,T}=nothing)
# Cache-less but non-allocating if `fx` and `epsilon` are supplied
# fx must be f(x)
DiffEqDiffTools.finite_difference_derivative(
f,
x :: AbstractArray{<:Number},
fdtype :: Type{T1} = Val{:central},
returntype :: Type{T2} = eltype(x), # return type of f
fx :: Union{Nothing,AbstractArray{<:Number}} = nothing,
epsilon :: Union{Nothing,AbstractArray{<:Real}} = nothing;
[epsilon_factor])
DiffEqDiffTools.finite_difference_derivative!(
df :: AbstractArray{<:Number},
f,
x :: AbstractArray{<:Number},
fdtype :: Type{T1} = Val{:central},
returntype :: Type{T2} = eltype(x),
fx :: Union{Nothing,AbstractArray{<:Number}} = nothing,
epsilon :: Union{Nothing,AbstractArray{<:Real}} = nothing;
[epsilon_factor])
# Cached
DiffEqDiffTools.finite_difference_derivative!(
df::AbstractArray{<:Number},
f,
x::AbstractArray{<:Number},
cache::DerivativeCache{T1,T2,fdtype,returntype};
[epsilon_factor])
DiffEqDiffTools.DerivativeCache(
x :: AbstractArray{<:Number},
fx :: Union{Nothing,AbstractArray{<:Number}} = nothing,
epsilon :: Union{Nothing,AbstractArray{<:Real}} = nothing,
fdtype :: Type{T1} = Val{:central},
returntype :: Type{T2} = eltype(x))
This allocates either fx
or epsilon
if these are nothing and they are needed.
fx
is the current call of f(x)
and is required for forward-differencing
(otherwise is not necessary).
Gradients are either a vector->scalar map f(x)
, or a scalar->vector map
f(fx,x)
if inplace=Val{true}
and fx=f(x)
if inplace=Val{false}
.
# Cache-less
DiffEqDiffTools.finite_difference_gradient(
f,
x,
fdtype::Type{T1}=Val{:central},
returntype::Type{T2}=eltype(x),
inplace::Type{Val{T3}}=Val{true};
[epsilon_factor])
DiffEqDiffTools.finite_difference_gradient!(
df,
f,
x,
fdtype::Type{T1}=Val{:central},
returntype::Type{T2}=eltype(df),
inplace::Type{Val{T3}}=Val{true};
[epsilon_factor])
# Cached
DiffEqDiffTools.finite_difference_gradient!(
df::AbstractArray{<:Number},
f,
x::AbstractArray{<:Number},
cache::GradientCache;
[epsilon_factor])
DiffEqDiffTools.GradientCache(
df :: Union{<:Number,AbstractArray{<:Number}},
x :: Union{<:Number, AbstractArray{<:Number}},
fdtype :: Type{T1} = Val{:central},
returntype :: Type{T2} = eltype(df),
inplace :: Type{Val{T3}} = Val{true})
DiffEqDiffTools.GradientCache(
c1 :: Union{Nothing,AbstractArray{<:Number}},
c2 :: Union{Nothing,AbstractArray{<:Number}},
fx :: Union{Nothing,<:Number,AbstractArray{<:Number}} = nothing,
fdtype :: Type{T1} = Val{:central},
returntype :: Type{T2} = eltype(df),
inplace :: Type{Val{T3}} = Val{true})
Note that here fx
is a cached function call of f
. If you provide fx
, then
fx
will be used in the forward differencing method to skip a function call.
It is on you to make sure that you update cache.fx
every time before
calling DiffEqDiffTools.finite_difference_gradient!
. A good use of this is if you have a
cache array for the output of fx
already being used, you can make it alias
into the differencing algorithm here.
Jacobians are for functions f!(fx,x)
when using in-place finite_difference_jacobian!
,
and fx = f(x)
when using out-of-place finite_difference_jacobain
. The out-of-place
jacobian will return a similar type as jac_prototype
if it is not a nothing
. For non-square
Jacobians, a cache which specifies the vector fx
is required.
For sparse differentiation, pass a colorvec
of matrix colors. sparsity
should be a sparse
or structured matrix (Tridiagonal
, Banded
, etc. according to the ArrayInterface.jl specs)
to allow for decompression, otherwise the result will be the colorvec compressed Jacobian.
# Cache-less
DiffEqDiffTools.finite_difference_jacobian(
f,
x :: AbstractArray{<:Number},
fdtype :: Type{T1}=Val{:central},
returntype :: Type{T2}=eltype(x),
relstep=default_relstep(fdtype, eltype(x)),
absstep=relstep,
colorvec = 1:length(x),
sparsity = nothing,
jac_prototype = nothing)
finite_difference_jacobian!(J::AbstractMatrix,
f,
x::AbstractArray{<:Number},
fdtype :: Type{T1}=Val{:forward},
returntype :: Type{T2}=eltype(x),
f_in :: Union{T2,Nothing}=nothing;
relstep=default_relstep(fdtype, eltype(x)),
absstep=relstep,
colorvec = 1:length(x),
sparsity = ArrayInterface.has_sparsestruct(J) ? J : nothing)
# Cached
DiffEqDiffTools.finite_difference_jacobian(
f,
x,
cache::JacobianCache;
relstep=default_relstep(fdtype, eltype(x)),
absstep=relstep,
colorvec = cache.colorvec,
sparsity = cache.sparsity,
jac_prototype = nothing)
DiffEqDiffTools.finite_difference_jacobian!(
J::AbstractMatrix{<:Number},
f,
x::AbstractArray{<:Number},
cache::JacobianCache;
relstep=default_relstep(fdtype, eltype(x)),
absstep=relstep,
colorvec = cache.colorvec,
sparsity = cache.sparsity)
DiffEqDiffTools.JacobianCache(
x,
fdtype :: Type{T1} = Val{:central},
returntype :: Type{T2} = eltype(x),
colorvec = 1:length(x)
sparsity = nothing)
This assumes the Jacobian is square.
DiffEqDiffTools.JacobianCache(
x1 ,
fx ,
fx1,
fdtype :: Type{T1} = Val{:central},
returntype :: Type{T2} = eltype(fx),
colorvec = 1:length(x),
sparsity = nothing)
Hessians are for functions f(x)
which return a scalar.
#Cacheless
finite_difference_hessian(f, x::AbstractArray{<:Number},
fdtype :: Type{T1}=Val{:hcentral},
inplace :: Type{Val{T2}} = x isa StaticArray ? Val{true} : Val{false};
relstep=default_relstep(fdtype, eltype(x)),
absstep=relstep)
finite_difference_hessian!(H::AbstractMatrix,f,
x::AbstractArray{<:Number},
fdtype :: Type{T1}=Val{:hcentral},
inplace :: Type{Val{T2}} = x isa StaticArray ? Val{true} : Val{false};
relstep=default_relstep(fdtype, eltype(x)),
absstep=relstep)
#Cached
finite_difference_hessian(
f,x,
cache::HessianCache{T,fdtype,inplace};
relstep=default_relstep(fdtype, eltype(x)),
absstep=relstep)
finite_difference_hessian!(H,f,x,
cache::HessianCache{T,fdtype,inplace};
relstep = default_relstep(fdtype, eltype(x)),
absstep = relstep)
HessianCache(x,fdtype::Type{T1}=Val{:hcentral},
inplace::Type{Val{T2}} = x isa StaticArray ? Val{true} : Val{false})
HessianCache(xpp,xpm,xmp,xmm,
fdtype::Type{T1}=Val{:hcentral},
inplace::Type{Val{T2}} = x isa StaticArray ? Val{true} : Val{false})
01/11/2017
9 days ago
204 commits