# Best practices for automatic generation of Jacobians

**URL:** <https://discourse.julialang.org/t/best-practices-for-automatic-generation-of-jacobians/68035>\
**Category:** General Usage\
**Tags:** pde, forwarddiff\
**Created:** [September 11, 2021, 8:47pm UTC](https://discourse.julialang.org/t/best-practices-for-automatic-generation-of-jacobians/68035 "2021-09-11T20:47:56Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Gravlax](https://avatars.discourse-cdn.com/v4/letter/g/dbc845/32.png) [@Gravlax](https://discourse.julialang.org/u/Gravlax)\
**Post date:** [September 11, 2021, 8:47pm UTC](https://discourse.julialang.org/t/best-practices-for-automatic-generation-of-jacobians/68035/1 "2021-09-11T20:47:56Z")

</div>

Hi everyone,

I would like to discuss the best practices for automatic generation of Jacobians when solving PDEs. I understand that the ‘way to go’ is to (1) define a residual function (in loop-style - pure Julia), (2) use SparsityDetection to get the pattern (assuming we don’t know a priori), (3) generate coloring and (4) finally use forwarddiff.  
Let’s say I want to solve a 1D Poisson with variable coefficient. So, I would start by defining a residual function. For automatic differentiation routines the number of input parameters is somehow limited. I manage to get it working with either defining parameters (like the coefficient and grid spacing) as global variables which is not a very satisfying approach (see MWE1, just below).

```julia
using ForwardDiff, SparseDiffTools, SparseArrays, LinearAlgebra, SparsityDetection
# These parameters that should not not be defined globally ideally
nx = 11
k = ones(Float64,nx+1) 
dx = 1.0
T = ones(Float64,nx) 
##################################################
function Func!( f, T ) # 2 arguments needed, others can be made optional
    for i in 2:length(T)-1 # a loop seems essential here, the vectorised version does not appear to work
        qW = -k[i] * (T[i] - T[i-1]) / dx
        qE = -k[i+1] * (T[i+1] - T[i]) / dx
        f[i] = (qE - qW) / dx
    end
    f[1] = -2T[1] + T[2]
    f[end] = T[end-1] - 2T[end]
    nothing
end
##################################################
function main()  
    # Sparsity pattern
    input = rand(nx)
    output = similar(input)
    sparsity = jacobian_sparsity(Func!,output,input)
    jac = Float64.(sparse(sparsity))
    # Makes coloring
    colors = matrix_colors(jac) 
    # Evaluate Jacobian
    forwarddiff_color_jacobian!(jac, Func!, T, colorvec = colors)
    # Compute residual
    f = zero(T)
    Func!( f, T )
    println("Before solve: norm f = ", norm(f))
    # Newton step
    dT = -jac\f
    T .+= dT
    # Compute residual
    Func!( f, T )
    println("After solve: norm f = ", norm(f))
end
##################################################
main()

```

Another solution is to define all input parameters in the main routine (what I’d like to do), then these parameters need to be parsed in my residual evaluation routine. In order to work, these parameters need to be assigned to default dummy values which sounds a bit awkward. (these are not useful for the differentiation itself as long as they do not depend on the solution, but are essential for direct function evaluations). This is what is done in MWE2 below.

```julia
using ForwardDiff, SparseDiffTools, SparseArrays, LinearAlgebra, SparsityDetection
##################################################
function Func!( f, T, k=ones(Float64, length(T)+1), dx=1.0 ) # 2 arguments needed, others can be made optional
    for i in 2:length(T)-1 # a loop seems essential here, the vectorised version does not appear to work
        qW = -k[i] * (T[i] - T[i-1]) / dx
        qE = -k[i+1] * (T[i+1] - T[i]) / dx
        f[i] = (qE - qW) / dx
    end
    f[1] = -2T[1] + T[2]
    f[end] = T[end-1] - 2T[end]
    nothing
end
##################################################
function main()
    # Configuration
    # These parameters that should not not be defined globally ideally
    nx = 11
    k = ones(Float64,nx+1) 
    dx = 1.0
    T = ones(Float64,nx)   
    # Sparsity pattern
    input = rand(nx)
    output = similar(input)
    sparsity = jacobian_sparsity(Func!,output,input)
    jac = Float64.(sparse(sparsity))
    # Makes coloring
    colors = matrix_colors(jac) 
    # Evaluate Jacobian
    forwarddiff_color_jacobian!(jac, Func!, T, colorvec = colors)
    # Compute residual
    f = zero(T)
    Func!( f, T, k, dx )
    println("Before solve: norm f = ", norm(f))
    # Newton step
    dT = -jac\f
    T .+= dT
    # Compute residual
    Func!( f, T, k, dx )
    println("After solve: norm f = ", norm(f))
end
##################################################
main()

```

So the question is, are there better ways to define the residual functions and evaluate the Jacobian? I would like to know are the best practices to then further move to multi-dimensional problems including non-linearities…

I realize that this topic is close to [that one](https://discourse.julialang.org/t/forwarddiff-complicated-functions/31094/13) but I’ve somehow not been able to apply the tricks discussed there to my example…

Thanks in advance for your hints!  
Cheers

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 11, 2021, 10:26pm UTC](https://discourse.julialang.org/t/best-practices-for-automatic-generation-of-jacobians/68035/2 "2021-09-11T22:26:54Z")

</div>

`jacobian_sparsity` splats the args, so you can just do `jacobian_sparsity(Func!,output,input, T, k, dx)` like the docs say and you should be fine?

---

<div class="post-metadata">

**Author:** ![Gravlax](https://avatars.discourse-cdn.com/v4/letter/g/dbc845/32.png) [@Gravlax](https://discourse.julialang.org/u/Gravlax)\
**Post date:** [September 12, 2021, 11:42am UTC](https://discourse.julialang.org/t/best-practices-for-automatic-generation-of-jacobians/68035/3 "2021-09-12T11:42:42Z")

</div>

Yes, thanks a lot! `jacobian_sparsity` works indeed well with splat arguments.  
Now the challenge is to get it to work with `forwarddiff_color_jacobian!`.  
So far it only works if ‘secondary’ arguments receive a default value like:  
`Func!( f, T, k=ones(length(T)+1), dx=1.0 )`  
This is maybe not a bad thing (?). I’ll upload an updated MWE if I find a better solution.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 12, 2021, 11:51am UTC](https://discourse.julialang.org/t/best-practices-for-automatic-generation-of-jacobians/68035/4 "2021-09-12T11:51:04Z")

</div>

Just use a closure. `(du,u) -> f(du,u,T,k,dx)`

---

<div class="post-metadata">

**Author:** ![Gravlax](https://avatars.discourse-cdn.com/v4/letter/g/dbc845/32.png) [@Gravlax](https://discourse.julialang.org/u/Gravlax)\
**Post date:** [September 15, 2021, 4:54pm UTC](https://discourse.julialang.org/t/best-practices-for-automatic-generation-of-jacobians/68035/5 "2021-09-15T16:54:31Z")

</div>

> [@ChrisRackauckas](#):
>
> (du,u) → f(du,u,T,k,dx)

Thank you!

The updated MWE:

```julia
using ForwardDiff, SparseDiffTools, SparseArrays, LinearAlgebra, SparsityDetection
##################################################
function Func!( f, T, k, dx ) # 2 arguments needed, others can be made optional
    for i in 2:length(T)-1 # a loop seems essential here, the vectorised version does not appear to work
        qW = -k[i] * (T[i] - T[i-1]) / dx
        qE = -k[i+1] * (T[i+1] - T[i]) / dx
        f[i] = (qE - qW) / dx
    end
    f[1] = -2T[1] + T[2] # Dirichlet (half dx away from BC)
    f[end] = T[end-1] - 2T[end] # Dirichlet
    nothing
end
##################################################
function main()
    # Configuration
    nx = 11
    k = ones(Float64,nx+1) 
    dx = 1.0
    T = ones(Float64,nx)  
    f = zero(T) 
    # Sparsity pattern
    input = rand(nx)
    output = similar(input)
    sparsity = jacobian_sparsity(Func!,output,input, k, dx)
    jac = Float64.(sparse(sparsity))
    # Makes coloring
    colors = matrix_colors(jac) 
    # Closure 
    f_closed = (f,T) -> Func!(f, T, k, dx)
    # Evaluate Jacobian
    forwarddiff_color_jacobian!(jac, f_closed, T, colorvec = colors)
    # Compute residual
    Func!( f, T, k, dx )
    println("Before solve: norm f = ", norm(f))
    # Newton step
    dT = -jac\f
    T .+= dT
    # Compute residual
    Func!( f, T, k, dx )
    println("After solve: norm f = ", norm(f))
end
##################################################
main()

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 27, 2022, 9:19pm UTC](https://discourse.julialang.org/t/best-practices-for-automatic-generation-of-jacobians/68035/6 "2022-08-27T21:19:26Z")

</div>

Does that mean it’s fine now?

---

<div class="post-metadata">

**Author:** ![Gravlax](https://avatars.discourse-cdn.com/v4/letter/g/dbc845/32.png) [@Gravlax](https://discourse.julialang.org/u/Gravlax)\
**Post date:** [August 27, 2022, 9:26pm UTC](https://discourse.julialang.org/t/best-practices-for-automatic-generation-of-jacobians/68035/7 "2022-08-27T21:26:48Z")

</div>

Well, it worked fine when I uploaded it a while ago. I realise now that I forgot to check the solution box. Not sure it is still up to date with latest development of the involved packages though… thanks again

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 27, 2022, 9:33pm UTC](https://discourse.julialang.org/t/best-practices-for-automatic-generation-of-jacobians/68035/8 "2022-08-27T21:33:24Z")

</div>

Okay cool, I wasn’t sure what this thread was about anymore 😅