Does NonlinearSolve use the analytical jacobian?

This question came up when I did some tests related to this post.

Imagine I have the code

using NonlinearSolve

function f(u, p)
    return u^2 - 2

function df(u, p)
    return 2u

fn = NonlinearFunction(f, jac=df)
prob = NonlinearProblem(fn, 1.0)
sol = solve(prob, NewtonRaphson(; concrete_jac=true))

Running this code does not print "Hey" to the terminal, so I guess that automatic differentiation is being used to compute the jacobian instead of the given function. I noticed this behavior because in the post that I mentioned, there was an incompatible interaction between AD and QuadGK that was present even if you provide the analytical jacobian.

Is this intended behavior? How could I provide the analytical jacobian to a problem? Also, does it need to be given to NonlinearFunction always and cannot be an optional parameter to NonlinearProblem?

At least in a quick review of the documentation I did not see a clear explanation on how to provide analytical jacobians.

PD: I am using the version

  [8913a72c] NonlinearSolve v3.13.1


I am experiencing the same issue. In my case, NonlinearSolve is used to find the root of a scalar equation f(u,p)=0. This equation appears in the evaluation of a constraint of an optimization problem that I am trying to solve with Optimization.jl.

A workaround seems to be rewriting the scalar function as a vector function, with u vector of length one. For example:

using NonlinearSolve

f_scl(u,p) = cos(u) - p*u
f_vec(u,p) = [cos(u[1]) - p*u[1]]

function J_scl(u,p)
    return -sin(u) - p

function J_vec(u,p)
    return [-sin(u[1]) - p;;]

u0 = 0.0
p = 3.0
alg = SimpleNewtonRaphson()

nf_scl = NonlinearFunction(f_scl,jac=J_scl)
pb_scl = NonlinearProblem(nf_scl,u0,p)
sl_scl = solve(pb_scl,alg)
println("scalar solution: ", sl_scl.u)

nf_vec = NonlinearFunction(f_vec,jac=J_vec)
pb_vec = NonlinearProblem(nf_vec,[u0],p)
sl_vec = solve(pb_vec,alg)
println("vector solution: ", sl_vec.u[1])

output of the above:

scalar solution: 0.3167508287712211
vector solution: 0.3167508287712211

so indeed the analytical Jacobian is evaluated only in the second case.
I tried to swap SimpleNewtonRaphson() with NewtonRaphson() but without luck.

Is there any way to force the solver to use the analytical Jacobian also in the scalar case?

Oh, so it is only in the scalar case, thank you for detecting that. I opened an issue on GitHub since this seems mostly a bug, because this behavior is most surely unintended.

for scalar functions if it is using a bisection type method (i.e. ITP), that would make a lot of sense. You don’t need jacobians for scalar equations.

Yeah, but what if you use NewtonRaphson? I would expect that if you provide the analytical derivative, then it uses it. Otherwise the interface is a little bit misleading.

Can you open an issue?

already exists: Analytical Jacobian is not being used in the 3.13 version of NonlinearSolve · Issue #451 · SciML/NonlinearSolve.jl · GitHub

1 Like