Pass Jacobian function to NLsolve

[briochemc]

I'm not sure but I don't know how to pass a Jacobian function to NLsolve when I solve via

solve(prob, SSRootfind(...))

It seems only the "state" function f! is grabbed from prob (and not the Jacobian, even if it is there). Is that correct? Any way to fix that so that one can then solve for steady-states using efficient Newton-like methods?

[ChrisRackauckas]

Yeah we'd need to modify https://github.com/SciML/SteadyStateDiffEq.jl/blob/master/src/solve.jl#L30 to handle this.

[briochemc]

Yes that's what I thought. I can try a PR :)

Side-ish question: Is there a reason why this __solve method dispatches on the abstract type alg::SteadyStateDiffEqAlgorithm and then checks if typeof(alg) <: SSRootfind?

Is that better than dispatching directly on alg::SSRootfind and having a fallback on the abstract type, i.e., something like __solve(..., alg::SteadyStateDiffEqAlgorithm, ...) = error("Algorithm not recognized")?

[ChrisRackauckas]

I thought it would evolve to support multiple libraries. Instead, in 2020 it should just use NonlinearProblem and take in a solver for that.

[briochemc]

I am confused 😅. I was talking about the dispatch on the algorithm, but you are mentioning the problem type. Does that mean that the AbstractSteadyStateProblems are going to be deprecated and you want that replaced by a AbstractNonLinearProblem?

[ChrisRackauckas]

No, I want SSRootfind to be a method that builds a NonlinearProblem and solves it, which would then work with all packages that can solve a NonlinearProblem.

[briochemc]

Ok I think I see. I will try to write something up taking this into account.

[briochemc]

Actually I'm still a bit confused. Sorry if I need much hand-holding, I'm trying! I'll split this post in two, API and Jacobian use, for clarity.

API

AFAIU, SSRootfind does not build a problem but just wraps a function that solves a problem given f, u0, and abstol: https://github.com/SciML/SteadyStateDiffEq.jl/blob/a841aecd577af55fe48806ec6468f9d571c6df0c/src/algorithms.jl#L6

So maybe I can try to outline an API and you can tell me where I'm confusing things 😅 I think I want to do something like

prob = SteadyStateProblem(f, u0, p) # where `f` is an `ODEfunction that may contain a jacobian`
sol = solve(prob, SSRootfind(alg))

My understanding is that you want this to be recast into

prob = NonlinearProblem(f, u0, p)
sol = solve(prob, alg)

So solve(prob, alg::SSRootfind) should be recast the SteadyStateProblem as a NonlinearProblem and do something like

function solve(prob, alg::SSRootfind)
    NLprob = NonlinearProblem(prob.f, prob.u0, prob.p)
    return solve(NLprob, alg.nlsolve)
end

Is that correct?

Use of the Jacobian

I guess that in place of https://github.com/SciML/SteadyStateDiffEq.jl/blob/a841aecd577af55fe48806ec6468f9d571c6df0c/src/solve.jl#L30 we could do something like (omitting the ! for readbility)

df = isnothing(f.jac) ? f : OnceDifferentiable(f.f, f.jac, ...)
u = alg.nlsolve(df, u0, abstol)

to let nlsolve know about the jacobian if it is provided, but I'm not sure... What do you think?

[ChrisRackauckas]

Is that correct?

Yup exactly.

we could do something like (omitting the ! for readbility)

We should expand the whole NonlinearProblem to be like ODEProblem: have a NonlinearFunction with a spot for the Jacobian and all of that. Then in this algorithm you just set all of those pieces. We'll want to at least have dispatches for nonlinear solving to Sundials (Kinsol) and NLsolve.jl, and other packages could follow after. This is a bit more work of course, but getting this part of the interface complete and inline with the differential equations will have some pretty nice advantages in the near future.

[ChrisRackauckas]

Deprecated for NonlinearSolve.jl's SciMLNLSolve