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

Specific Domains Optimization (Mathematical)
1

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);
2

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

1 Like
3

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
4

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)

5

Thank you for your answer,

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,

6

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
7

Thanks!