# Correct use of the sparsity pattern in NLsolve.jl for Diagonal Jacobians

Hello,

I would like to solve \sim 400 nonlinear equations using NLsolve.jl for which the Jacobian matrix is diagonal, i.e. the components are independent of each other.
I would like to use the sparsity pattern to speed up the resolution. I have looked at the file sparse.jl of the tests folder to construct my example. The function used is just a toy problem, in my real problem it’s a more complex function.

How can I optimize this code? Are there methods that support hessian as well?
Is it better to use the tools from SparsityArrays or to provide the Jacobian as a `Diagonal` object?

``````using SparseArrays
using NLsolve

Ne = 400

function f_sparse!(F, x)
@inbounds for (i, xi) in enumerate(x)
F[i] = exp(-xi^2/(2*i))*erf(xi)
end
end

function g_sparse!(J, x)
@inbounds for (i, xi) in enumerate(x)
J[i,i] = (-xi/i)*exp(-xi^2/(2*i))*erf(xi) + exp(-xi^2/(2*i))*(2/√π)*exp(-xi^2)
end
end

df = OnceDifferentiable(f_sparse!, g_sparse!, randn(Ne), rand(Ne), spzeros(Ne, Ne))

F = rand(Ne);

r = nlsolve(df, F, method = :newton);
``````

Provide the Jacobian as diagonal and importantly use SparseDiffTools to speed up the Jacobian construction.

1 Like

Thank you for your help. Do you know where I can find an example that use NLsolve and SparseDiffTools?

I provide now the Jacobian as a diagonal matrix but I don’t know how to use the package SparseDiffTools.jl

``````using NLsolve
using SparseArrays
using SparseDiffTools
using SparsityDetection

Ne = 200
function f!(F, x)
@inbounds for (i, xi) in enumerate(x)
F[i] = exp(-xi^2/(2*i))*erf(xi)
end
nothing
end

function g_diag!(J, x)
@inbounds for (i, xi) in enumerate(x)
J[i,i] = (-xi/i)*exp(-xi^2/(2*i))*erf(xi) + exp(-xi^2/(2*i))*(2/√π)*exp(-xi^2)
end
J = Diagonal(J)
nothing
end

function g_sparse!(J, x)
@inbounds for (i, xi) in enumerate(x)
J[i,i] = (-xi/i)*exp(-xi^2/(2*i))*erf(xi) + exp(-xi^2/(2*i))*(2/√π)*exp(-xi^2)
end
nothing
end

df_diag = OnceDifferentiable(f!, g_diag!, randn(Ne), rand(Ne), Diagonal(zeros(Ne)))
df_sparse = OnceDifferentiable(f!, g_sparse!, randn(Ne), rand(Ne), spzeros(Ne,Ne))

F = rand(Ne);

@btime r_diag = nlsolve(\$df_diag, \$F, method = :newton);
@btime r_sparse = nlsolve(\$df_sparse, \$F, method = :newton);
``````

124.516 μs (91 allocations: 95.00 KiB)
237.630 μs (511 allocations: 312.55 KiB)

I don’t know how to set the sparsity pattern to `Diagonal`, I got the following error message

``````input = zeros(Ne)
output = similar(input)
sparsity_pattern = jacobian_sparsity(f!,output,input)
jac = Float64.(sparse(sparsity_pattern));
``````
``````ArgumentError: invalid index: Tagged(Tag{nametype(JacobianSparsityContext),1031646979481733957,Nothing}(), 2, _) of type Cassette.Tagged{Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Int64,SparsityDetection.ProvinanceSet,Cassette.NoMetaMeta,Cassette.Context{nametype(JacobianSparsityContext),Tuple{Sparsity,SparsityDetection.Path},Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},SparsityDetection.var"##PassType#253",IdDict{Module,Dict{Symbol,Cassette.BindingMeta}},Cassette.DisableHooks}}

Stacktrace:
[1] to_index(::Cassette.Tagged{Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Int64,SparsityDetection.ProvinanceSet,Cassette.NoMetaMeta,Cassette.Context{nametype(JacobianSparsityContext),Tuple{Sparsity,SparsityDetection.Path},Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},SparsityDetection.var"##PassType#253",IdDict{Module,Dict{Symbol,Cassette.BindingMeta}},Cassette.DisableHooks}}) at ./indices.jl:297
[2] to_index(::LinearIndices{1,Tuple{Base.OneTo{Int64}}}, ::Cassette.Tagged{Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Int64,SparsityDetection.ProvinanceSet,Cassette.NoMetaMeta,Cassette.Context{nametype(JacobianSparsityContext),Tuple{Sparsity,SparsityDetection.Path},Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},SparsityDetection.var"##PassType#253",IdDict{Module,Dict{Symbol,Cassette.BindingMeta}},Cassette.DisableHooks}}) at ./indices.jl:274
[3] to_indices at ./indices.jl:325 [inlined]
[4] to_indices at ./indices.jl:322 [inlined]
[5] getindex at ./abstractarray.jl:980 [inlined]
[6] overdub at /home/mat/.julia/packages/SparsityDetection/PbKzm/src/jacobian.jl:98 [inlined]
[7] iterate at ./array.jl:763 [inlined]
[8] recurse(::Cassette.Context{nametype(JacobianSparsityContext),Tuple{Sparsity,SparsityDetection.Path},Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},SparsityDetection.var"##PassType#253",IdDict{Module,Dict{Symbol,Cassette.BindingMeta}},Cassette.DisableHooks}, ::typeof(iterate), ::Cassette.Tagged{Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Array{Float64,1},Union{SparsityDetection.JacInput, SparsityDetection.JacOutput, SparsityDetection.ProvinanceSet},Array{Cassette.Meta{SparsityDetection.ProvinanceSet,Cassette.NoMetaMeta},1},Cassette.Context{nametype(JacobianSparsityContext),Tuple{Sparsity,SparsityDetection.Path},Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},SparsityDetection.var"##PassType#253",IdDict{Module,Dict{Symbol,Cassette.BindingMeta}},Cassette.DisableHooks}}, ::Cassette.Tagged{Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Int64,SparsityDetection.ProvinanceSet,Cassette.NoMetaMeta,Cassette.Context{nametype(JacobianSparsityContext),Tuple{Sparsity,SparsityDetection.Path},Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},SparsityDetection.var"##PassType#253",IdDict{Module,Dict{Symbol,Cassette.BindingMeta}},Cassette.DisableHooks}}) at /home/mat/.julia/packages/Cassette/158rp/src/overdub.jl:0
[9] overdub(::Cassette.Context{nametype(JacobianSparsityContext),Tuple{Sparsity,SparsityDetection.Path},Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},SparsityDetection.var"##PassType#253",IdDict{Module,Dict{Symbol,Cassette.BindingMeta}},Cassette.DisableHooks}, ::Function, ::Cassette.Tagged{Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Array{Float64,1},Union{SparsityDetection.JacInput, SparsityDetection.JacOutput, SparsityDetection.ProvinanceSet},Array{Cassette.Meta{SparsityDetection.ProvinanceSet,Cassette.NoMetaMeta},1},Cassette.Context{nametype(JacobianSparsityContext),Tuple{Sparsity,SparsityDetection.Path},Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},SparsityDetection.var"##PassType#253",IdDict{Module,Dict{Symbol,Cassette.BindingMeta}},Cassette.DisableHooks}}, ::Cassette.Tagged{Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Int64,SparsityDetection.ProvinanceSet,Cassette.NoMetaMeta,Cassette.Context{nametype(JacobianSparsityContext),Tuple{Sparsity,SparsityDetection.Path},Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},SparsityDetection.var"##PassType#253",IdDict{Module,Dict{Symbol,Cassette.BindingMeta}},Cassette.DisableHooks}}) at /home/mat/.julia/packages/SparsityDetection/PbKzm/src/propagate_tags.jl:45
[10] overdub(::Cassette.Context{nametype(JacobianSparsityContext),Tuple{Sparsity,SparsityDetection.Path},Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},SparsityDetection.var"##PassType#253",IdDict{Module,Dict{Symbol,Cassette.BindingMeta}},Cassette.DisableHooks}, ::typeof(Core._apply_iterate), ::Function, ::Function, ::Cassette.Tagged{Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Tuple{Array{Float64,1}},Union{SparsityDetection.JacInput, SparsityDetection.JacOutput, SparsityDetection.ProvinanceSet},Tuple{Cassette.Immutable{Cassette.Meta{Union{SparsityDetection.JacInput, SparsityDetection.JacOutput, SparsityDetection.ProvinanceSet},Array{Cassette.Meta{SparsityDetection.ProvinanceSet,Cassette.NoMetaMeta},1}}}},Cassette.Context{nametype(JacobianSparsityContext),Tuple{Sparsity,SparsityDetection.Path},Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},SparsityDetection.var"##PassType#253",IdDict{Module,Dict{Symbol,Cassette.BindingMeta}},Cassette.DisableHooks}}, ::Cassette.Tagged{Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Tuple{Int64},Union{SparsityDetection.JacInput, SparsityDetection.JacOutput, SparsityDetection.ProvinanceSet},Tuple{Cassette.Immutable{Cassette.Meta{SparsityDetection.ProvinanceSet,Cassette.NoMetaMeta}}},Cassette.Context{nametype(JacobianSparsityContext),Tuple{Sparsity,SparsityDetection.Path},Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},SparsityDetection.var"##PassType#253",IdDict{Module,Dict{Symbol,Cassette.BindingMeta}},Cassette.DisableHooks}}) at /home/mat/.julia/packages/Cassette/158rp/src/context.jl:266
[11] iterate at ./iterators.jl:139 [inlined]
[12] recurse(::Cassette.Context{nametype(JacobianSparsityContext),Tuple{Sparsity,SparsityDetection.Path},Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},SparsityDetection.var"##PassType#253",IdDict{Module,Dict{Symbol,Cassette.BindingMeta}},Cassette.DisableHooks}, ::typeof(iterate), ::Cassette.Tagged{Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Base.Iterators.Enumerate{Array{Float64,1}},Union{SparsityDetection.JacInput, SparsityDetection.JacOutput, SparsityDetection.ProvinanceSet},NamedTuple{(:itr,),Tuple{Cassette.Immutable{Cassette.Meta{Union{SparsityDetection.JacInput, SparsityDetection.JacOutput, SparsityDetection.ProvinanceSet},Array{Cassette.Meta{SparsityDetection.ProvinanceSet,Cassette.NoMetaMeta},1}}}}},Cassette.Context{nametype(JacobianSparsityContext),Tuple{Sparsity,SparsityDetection.Path},Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},SparsityDetection.var"##PassType#253",IdDict{Module,Dict{Symbol,Cassette.BindingMeta}},Cassette.DisableHooks}}, ::Cassette.Tagged{Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Tuple{Int64,Int64},Union{SparsityDetection.JacInput, SparsityDetection.JacOutput, SparsityDetection.ProvinanceSet},Tuple{Cassette.Immutable{Cassette.Meta{SparsityDetection.ProvinanceSet,Cassette.NoMetaMeta}},Cassette.Immutable{Cassette.Meta{SparsityDetection.ProvinanceSet,Cassette.NoMetaMeta}}},Cassette.Context{nametype(JacobianSparsityContext),Tuple{Sparsity,SparsityDetection.Path},Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},SparsityDetection.var"##PassType#253",IdDict{Module,Dict{Symbol,Cassette.BindingMeta}},Cassette.DisableHooks}}) at /home/mat/.julia/packages/Cassette/158rp/src/overdub.jl:0
[13] overdub(::Cassette.Context{nametype(JacobianSparsityContext),Tuple{Sparsity,SparsityDetection.Path},Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},SparsityDetection.var"##PassType#253",IdDict{Module,Dict{Symbol,Cassette.BindingMeta}},Cassette.DisableHooks}, ::Function, ::Cassette.Tagged{Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Base.Iterators.Enumerate{Array{Float64,1}},Union{SparsityDetection.JacInput, SparsityDetection.JacOutput, SparsityDetection.ProvinanceSet},NamedTuple{(:itr,),Tuple{Cassette.Immutable{Cassette.Meta{Union{SparsityDetection.JacInput, SparsityDetection.JacOutput, SparsityDetection.ProvinanceSet},Array{Cassette.Meta{SparsityDetection.ProvinanceSet,Cassette.NoMetaMeta},1}}}}},Cassette.Context{nametype(JacobianSparsityContext),Tuple{Sparsity,SparsityDetection.Path},Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},SparsityDetection.var"##PassType#253",IdDict{Module,Dict{Symbol,Cassette.BindingMeta}},Cassette.DisableHooks}}, ::Cassette.Tagged{Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Tuple{Int64,Int64},Union{SparsityDetection.JacInput, SparsityDetection.JacOutput, SparsityDetection.ProvinanceSet},Tuple{Cassette.Immutable{Cassette.Meta{SparsityDetection.ProvinanceSet,Cassette.NoMetaMeta}},Cassette.Immutable{Cassette.Meta{SparsityDetection.ProvinanceSet,Cassette.NoMetaMeta}}},Cassette.Context{nametype(JacobianSparsityContext),Tuple{Sparsity,SparsityDetection.Path},Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},SparsityDetection.var"##PassType#253",IdDict{Module,Dict{Symbol,Cassette.BindingMeta}},Cassette.DisableHooks}}) at /home/mat/.julia/packages/SparsityDetection/PbKzm/src/propagate_tags.jl:37
[14] f! at ./In[113]:9 [inlined]
[15] recurse(::Cassette.Context{nametype(JacobianSparsityContext),Tuple{Sparsity,SparsityDetection.Path},Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},SparsityDetection.var"##PassType#253",IdDict{Module,Dict{Symbol,Cassette.BindingMeta}},Cassette.DisableHooks}, ::typeof(f!), ::Cassette.Tagged{Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Array{Float64,1},Union{SparsityDetection.JacInput, SparsityDetection.JacOutput, SparsityDetection.ProvinanceSet},Array{Cassette.Meta{SparsityDetection.ProvinanceSet,Cassette.NoMetaMeta},1},Cassette.Context{nametype(JacobianSparsityContext),Sparsity,Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Cassette.var"##PassType#255",IdDict{Module,Dict{Symbol,Cassette.BindingMeta}},Cassette.DisableHooks}}, ::Cassette.Tagged{Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Array{Float64,1},Union{SparsityDetection.JacInput, SparsityDetection.JacOutput, SparsityDetection.ProvinanceSet},Array{Cassette.Meta{SparsityDetection.ProvinanceSet,Cassette.NoMetaMeta},1},Cassette.Context{nametype(JacobianSparsityContext),Sparsity,Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Cassette.var"##PassType#255",IdDict{Module,Dict{Symbol,Cassette.BindingMeta}},Cassette.DisableHooks}}) at /home/mat/.julia/packages/Cassette/158rp/src/overdub.jl:0
[16] overdub(::Cassette.Context{nametype(JacobianSparsityContext),Tuple{Sparsity,SparsityDetection.Path},Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},SparsityDetection.var"##PassType#253",IdDict{Module,Dict{Symbol,Cassette.BindingMeta}},Cassette.DisableHooks}, ::Function, ::Cassette.Tagged{Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Array{Float64,1},Union{SparsityDetection.JacInput, SparsityDetection.JacOutput, SparsityDetection.ProvinanceSet},Array{Cassette.Meta{SparsityDetection.ProvinanceSet,Cassette.NoMetaMeta},1},Cassette.Context{nametype(JacobianSparsityContext),Sparsity,Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Cassette.var"##PassType#255",IdDict{Module,Dict{Symbol,Cassette.BindingMeta}},Cassette.DisableHooks}}, ::Cassette.Tagged{Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Array{Float64,1},Union{SparsityDetection.JacInput, SparsityDetection.JacOutput, SparsityDetection.ProvinanceSet},Array{Cassette.Meta{SparsityDetection.ProvinanceSet,Cassette.NoMetaMeta},1},Cassette.Context{nametype(JacobianSparsityContext),Sparsity,Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Cassette.var"##PassType#255",IdDict{Module,Dict{Symbol,Cassette.BindingMeta}},Cassette.DisableHooks}}) at /home/mat/.julia/packages/SparsityDetection/PbKzm/src/propagate_tags.jl:45
[17] overdub(::Cassette.Context{nametype(JacobianSparsityContext),Tuple{Sparsity,SparsityDetection.Path},Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},SparsityDetection.var"##PassType#253",IdDict{Module,Dict{Symbol,Cassette.BindingMeta}},Cassette.DisableHooks}, ::typeof(Core._apply_iterate), ::Function, ::Function, ::Tuple{Cassette.Tagged{Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Array{Float64,1},Union{SparsityDetection.JacInput, SparsityDetection.JacOutput, SparsityDetection.ProvinanceSet},Array{Cassette.Meta{SparsityDetection.ProvinanceSet,Cassette.NoMetaMeta},1},Cassette.Context{nametype(JacobianSparsityContext),Sparsity,Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Cassette.var"##PassType#255",IdDict{Module,Dict{Symbol,Cassette.BindingMeta}},Cassette.DisableHooks}},Cassette.Tagged{Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Array{Float64,1},Union{SparsityDetection.JacInput, SparsityDetection.JacOutput, SparsityDetection.ProvinanceSet},Array{Cassette.Meta{SparsityDetection.ProvinanceSet,Cassette.NoMetaMeta},1},Cassette.Context{nametype(JacobianSparsityContext),Sparsity,Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Cassette.var"##PassType#255",IdDict{Module,Dict{Symbol,Cassette.BindingMeta}},Cassette.DisableHooks}}}) at /home/mat/.julia/packages/Cassette/158rp/src/context.jl:266
[18] #2 at /home/mat/.julia/packages/SparsityDetection/PbKzm/src/controlflow.jl:148 [inlined]
[19] recurse(::Cassette.Context{nametype(JacobianSparsityContext),Tuple{Sparsity,SparsityDetection.Path},Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},SparsityDetection.var"##PassType#253",IdDict{Module,Dict{Symbol,Cassette.BindingMeta}},Cassette.DisableHooks}, ::SparsityDetection.var"#2#3"{typeof(f!),Tuple{Cassette.Tagged{Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Array{Float64,1},Union{SparsityDetection.JacInput, SparsityDetection.JacOutput, SparsityDetection.ProvinanceSet},Array{Cassette.Meta{SparsityDetection.ProvinanceSet,Cassette.NoMetaMeta},1},Cassette.Context{nametype(JacobianSparsityContext),Sparsity,Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Cassette.var"##PassType#255",IdDict{Module,Dict{Symbol,Cassette.BindingMeta}},Cassette.DisableHooks}},Cassette.Tagged{Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Array{Float64,1},Union{SparsityDetection.JacInput, SparsityDetection.JacOutput, SparsityDetection.ProvinanceSet},Array{Cassette.Meta{SparsityDetection.ProvinanceSet,Cassette.NoMetaMeta},1},Cassette.Context{nametype(JacobianSparsityContext),Sparsity,Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Cassette.var"##PassType#255",IdDict{Module,Dict{Symbol,Cassette.BindingMeta}},Cassette.DisableHooks}}}}) at /home/mat/.julia/packages/Cassette/158rp/src/overdub.jl:0
[20] abstract_run(::SparsityDetection.var"#22#24", ::Cassette.Context{nametype(JacobianSparsityContext),Sparsity,Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Cassette.var"##PassType#255",IdDict{Module,Dict{Symbol,Cassette.BindingMeta}},Cassette.DisableHooks}, ::Function, ::Cassette.Tagged{Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Array{Float64,1},Union{SparsityDetection.JacInput, SparsityDetection.JacOutput, SparsityDetection.ProvinanceSet},Array{Cassette.Meta{SparsityDetection.ProvinanceSet,Cassette.NoMetaMeta},1},Cassette.Context{nametype(JacobianSparsityContext),Sparsity,Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Cassette.var"##PassType#255",IdDict{Module,Dict{Symbol,Cassette.BindingMeta}},Cassette.DisableHooks}}, ::Vararg{Cassette.Tagged{Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Array{Float64,1},Union{SparsityDetection.JacInput, SparsityDetection.JacOutput, SparsityDetection.ProvinanceSet},Array{Cassette.Meta{SparsityDetection.ProvinanceSet,Cassette.NoMetaMeta},1},Cassette.Context{nametype(JacobianSparsityContext),Sparsity,Cassette.Tag{nametype(JacobianSparsityContext),0x0e512576419d0b45,Nothing},Cassette.var"##PassType#255",IdDict{Module,Dict{Symbol,Cassette.BindingMeta}},Cassette.DisableHooks}},N} where N; verbose::Bool) at /home/mat/.julia/packages/SparsityDetection/PbKzm/src/controlflow.jl:148
[21] jacobian_sparsity(::Function, ::Array{Float64,1}, ::Array{Float64,1}; sparsity::Sparsity, verbose::Bool, raw::Bool) at /home/mat/.julia/packages/SparsityDetection/PbKzm/src/jacobian.jl:135
[22] jacobian_sparsity(::Function, ::Array{Float64,1}, ::Array{Float64,1}) at /home/mat/.julia/packages/SparsityDetection/PbKzm/src/jacobian.jl:130
[23] top-level scope at In[115]:3
``````

You don’t have to use automated Jacobian sparsity discovery if you already know the pattern is Diagonal. The colorvec is just `ones(n)`, define the Jacobian to use FiniteDiff.jl or SparseDiffTools.jl with that colorvec and the right `sparsity`, and open an issue to let NLsolve let you choose the matrix type for the last optimization (lu)

I have modified my code to better use that the Jacobian is diagonal, in particular `g_diag!`. I am passing an extra argument `linsolve` to invert the diagonal Jacobian. I don’t understand how to use SparseDiffTools to speed-up Jacobian computations. I don’t need finite-difference or automatic differentiation, I can compute analytically the Jacobian.

I am not sure I understand, can you develop a bit

``````using NLsolve
using SparseArrays
using SparseDiffTools
using SparsityDetection

Ne = 200

function f!(F, x)
@inbounds for (i, xi) in enumerate(x)
F[i] = exp(-xi^2/(2*i))*erf(xi)
end
nothing
end

function g_diag!(J, x)
@inbounds for (i, xi) in enumerate(x)
J.diag[i] = (-xi/i)*exp(-xi^2/(2*i))*erf(xi) + exp(-xi^2/(2*i))*(2/√π)*exp(-xi^2)
end
nothing
end

function g_sparse!(J, x)
@inbounds for (i, xi) in enumerate(x)
J[i,i] = (-xi/i)*exp(-xi^2/(2*i))*erf(xi) + exp(-xi^2/(2*i))*(2/√π)*exp(-xi^2)
end
end

df_diag = OnceDifferentiable(f!, g_diag!, randn(Ne), rand(Ne), Diagonal(zeros(Ne)))
df_sparse = OnceDifferentiable(f!, g_sparse!, randn(Ne), rand(Ne), spzeros(Ne,Ne))

F = rand(Ne);

# Method with a Diagonal Jacobian
@btime r_diag = nlsolve(\$df_diag, \$F; linsolve = (x, A, b) -> copyto!(x, b ./ A.diag), method = :newton);

# Method with a sparse Jacobian
@btime r_sparse = nlsolve(\$df_sparse, \$F; method = :newton);
``````

114.696 μs (88 allocations: 91.44 KiB)
221.794 μs (488 allocations: 298.63 KiB)

``````input = zeros(Ne)
output = similar(input)
sparsity_pattern = ones(Ne)
jac = Float64.(sparse(sparsity_pattern))
``````

Thanks,

Oh, it’s for the case where the analytical Jacobian isn’t know. If you do know the analytical Jacobian, then you’re good.

1 Like

Thanks!