JuliaDiff / SparseDiffTools.jl

Fast jacobian computation through sparsity exploitation and matrix coloring
MIT License
237 stars 41 forks source link

Add method to mirror ForwardDiff's syntax for writing Jacobian of an oop function to an allocated matrix #128

Closed chenwilliam77 closed 3 years ago

chenwilliam77 commented 3 years ago

In ForwardDiff.jl, there is the ability to define an out-of-place function f(x) and autodiff it using

ForwardDiff.jacobian!(jac, f, x)

where jac is allocated outside of ForwardDiff.jacobian!. I have a specific use case where it'd be nice to have such a functionality in SparseDiffTools.jl and have implemented it. The tests should pass.

You may notice that I've defined

function forwarddiff_color_jacobian(J::SparseMatrixCSC{<:Number},f,x::AbstractArray{<:Number},jac_cache::ForwardColorJacCache,jac_prototype=nothing)

in addition to

function forwarddiff_color_jacobian(J::AbstractArray{<:Number},f,x::AbstractArray{<:Number},jac_cache::ForwardColorJacCache,jac_prototype=nothing)

I'm not too familiar with StaticArrays.jl, but I was unable to get one of the StaticArrays tests to pass when I used

J .+= (size(Ji)!=size(J) ? reshape(Ji,size(J)) : Ji)

instead of

J = J + (size(Ji)!=size(J) ? reshape(Ji,size(J)) : Ji)

To ensure I can mirror the function from ForwardDiff, I need the first line of code to work, but it seems to cause an error for StaticArrays.

Thus, to avoid changing the behavior of the package by too much, I just defined a version of forwarddiff_color_jacobian specifically when J is a sparse matrix, since it should be safe to do J .+= ... The more general version of forwarddiff_color_jacobian continues to use J += ... If you think that it'd be better to switch the two (i.e. write a forwarddiff_color_jacobian(J::StaticArray{<: Number}, ...) and using J .+= .. for the general version), then I'm happy to update the code.

codecov-io commented 3 years ago

Codecov Report

Merging #128 (9b24650) into master (9637980) will increase coverage by 1.10%. The diff coverage is 93.10%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master     #128      +/-   ##
==========================================
+ Coverage   80.46%   81.56%   +1.10%     
==========================================
  Files          12       12              
  Lines         558      613      +55     
==========================================
+ Hits          449      500      +51     
- Misses        109      113       +4     
Impacted Files Coverage Δ
src/differentiation/compute_jacobian_ad.jl 92.39% <93.10%> (+0.15%) :arrow_up:

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update 9637980...9b24650. Read the comment docs.

ChrisRackauckas commented 3 years ago

It might be good to separate the StaticArray case anyways by using ArrayInterface.ismutable(x). That way we can have a fully non-allocating version and a fully non-mutating version.

ChrisRackauckas commented 3 years ago

If SparseMatrixCSC has its own branch, then https://github.com/JuliaDiff/SparseDiffTools.jl/pull/128/files#diff-7e39b900cec80b7d52f89704e277c0f56afc3c2197d5ba4fef1fa9c943bcb0a2R135-R136 can be removed.

ChrisRackauckas commented 3 years ago

To ensure I can mirror the function from ForwardDiff, I need the first line of code to work, but it seems to cause an error for StaticArrays.

StaticArrays cannot be mutated, so they need a fully non-mutating path. We should simply throw an error if ArrayInterface.mutable(J) == false, since that means J[i] = ... will error.

ChrisRackauckas commented 3 years ago

Thus, to avoid changing the behavior of the package by too much, I just defined a version of forwarddiff_color_jacobian specifically when J is a sparse matrix, since it should be safe to do J .+= ... The more general version of forwarddiff_color_jacobian continues to use J += ... If you think that it'd be better to switch the two (i.e. write a forwarddiff_color_jacobian(J::StaticArray{<: Number}, ...) and using J .+= .. for the general version), then I'm happy to update the code.

We should do it slightly the other way around. Your version is a bit faster, so we should make your version be the one that is normally used, and we should send anything with ArrayInterface.mutable(J) == false to use the old OOP version since that would be slower for example on a dense J.

ChrisRackauckas commented 3 years ago

great work, thanks!