Does NonlinearSolve use the analytical jacobian?

General Usage
1

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
end

function df(u, p)
    println("Hey")
    return 2u
end

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
2

Hello,

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)
    println("hop!")
    return -sin(u) - p
end

function J_vec(u,p)
    println("hop!")
    return [-sin(u[1]) - p;;]
end

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
hop!
hop!
hop!
hop!
hop!
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?
Thanks!

3

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.

4

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.

5

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.

6

Can you open an issue?

7

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

1 Like