Error when using a neural network inside a system of nonlinear equations.

[GevikUCL]

I've set up the following very simple system of nonlinear equations to solve.

using NonlinearSolve, Lux, LinearAlgebra, StableRNGs

rng = StableRNG(1111)

function func(u, p)
    [u[1].*u[1] .- p[1],
    u[2].*u[2] .- p[2]]
end

guess = [1., 1.]
params = [2.,3.]

prob = NonlinearProblem{false}(func, guess, params)
solution = solve(prob, NewtonRaphson())

This works fine. I am now trying to learn one of the terms from this system of equations using a neural network (similar to the UDE approach https://docs.sciml.ai/Overview/dev/showcase/missing_physics/).

I define a network.

ann = Lux.Chain(
    Lux.Dense(2,20,tanh), Lux.Dense(20,20,tanh), Lux.Dense(20,1)
)
p, st = Lux.setup(rng, ann)

And define the grey-box model, where I replace the u[2]*u[2] term with the output of the network.

function GreyBox(u, p)
    z = ann(u, p, st)[1]
    [u[1].*u[1] .- params[1],
    z[1] .- params[2]]
end

NN_prob = NonlinearProblem{false}(GreyBox, guess, p)
NN_solution = solve(NN_prob NewtonRaphson())

In this case, solve isn't able to find the solution and the following error is returned.

ERROR: SingularException(2)

Was wondering if there would be a way to work around this issue.

[ChrisRackauckas]

What's the error message? Please share the whole message.

[ChrisRackauckas]

This seems to be because it is just singular at that point.

J = [2.0 0.0; 0.015492069072784054 -0.2379234976173196]
J = [3.0 0.0; 0.04881018347278711 0.011015270204509993]
J = [2.8333333333333335 0.0; 0.005072970655426244 -0.00010978688946150561]
J = [2.8284313725490198 0.0; 2.339615366002708e-18 -5.4331830373905e-20]
J = [2.8284271247493797 0.0; 0.0 0.0]

That's something that can come up with a Newton method. Meanwhile the following works:

NN_solution = solve(NN_prob, TrustRegion(5))

So this is a good test case for adding trust regioning to NewtonRaphson. Pinging @CCsimon123 so he's aware of this test case (and to note he's already solved some issues!).

[Deltadahl]

I will make a PR tomorrow with a new TrustRegion solver to NonlinearSolve.jl
It will hopefully handle this by default, but I will make sure to test it!

[ChrisRackauckas]

Indeed, this is just because a standard Newton cannot handle the singularity, but trust region can.

using NonlinearSolve, Lux, LinearAlgebra, StableRNGs

rng = StableRNG(1111)

function func(u, p)
    [u[1].*u[1] .- p[1],
    u[2].*u[2] .- p[2]]
end

guess = [1., 1.]
params = [2.,3.]

prob = NonlinearProblem{false}(func, guess, params)
solution = solve(prob)

ann = Lux.Chain(
    Lux.Dense(2,20,tanh), Lux.Dense(20,20,tanh), Lux.Dense(20,1)
)
p, st = Lux.setup(rng, ann)

function GreyBox(u, p)
    z = ann(u, p, st)[1]
    [u[1].*u[1] .- params[1],
    z[1] .- params[2]]
end

NN_prob = NonlinearProblem{false}(GreyBox, guess, p)
NN_solution = solve(NN_prob)

this is fine with the default methods we have now, closing.