Solving a nonlinear equation with NonlinearSolve.jl 2.X

Specific Domains Optimization (Mathematical)
1

Hello,

starting from NonlinearSolve.jl version 2.0 my code stopped working. I reduced it to the following MWE:

using NonlinearSolve

function exec_mwe()
    f(u, p) = [(u[1] - p[1])^2 + (u[2] - p[2])^2]
    u0 = rand(2)
    p = [1.0, 2.0]
    prob = NonlinearProblem(f, u0, p)
    sol = solve(prob, NewtonRaphson())
    println(sol)
end

exec_mwe()

This is the error message:

ERROR: DimensionMismatch: matrix is not square: dimensions are (1, 2)
Stacktrace:
[1] checksquare
@ C:\Users\np.julia\juliaup\julia-1.9.3+0.x64.w64.mingw32\share\julia\stdlib\v1.9\LinearAlgebra\src\LinearAlgebra.jl:239 [inlined]
[2] UpperTriangular
@ C:\Users\np.julia\juliaup\julia-1.9.3+0.x64.w64.mingw32\share\julia\stdlib\v1.9\LinearAlgebra\src\triangular.jl:17 [inlined]
…

I realize that there is a breaking change in NonlinearSolve.jl. However, I dont know how to adjust my code. Now it seems to be necessary to have the same number of equations and variables.

2

What system are you trying to solve?

In this line

you define a function that takes two inputs but only returns one output. In order to solve a system, you should have as many variables as equations, which is exactly what the error message is telling you.

1 Like
3

What you are saying is true for linear systems of equations. However, for nonlinear systems of equations it is not true in general. The example from the MWE has a unique solution u[1] = p[1] and u[2] = p[2] eventhough there are two variables and only one equation. In version 1.10.1 NonlinearSolve.jl was able to find this solution.

4

I think there may was something different in that formulation than what you are using right now. If you are using a Newton-Raphson algorithm, you use the linear approximation to your nonlinear function to approximate the real root from an initial guess. This jacobian should be a square matrix (which is not in your case) and then you are indeed solving linear systems repeteadly until you get convergence.

So if your nonlinear function has not as many variables as equations your jacobian is not square and that is why you get that error. You probably used another technique to solve that nonlinear equation in the past.

5

The MWE I posted works in 1.10.1 (with the Newton-Raphson-Method) and stopped working in 2.0. The condition that the Jacobian is a square matrix is new.

6

Define it as a NonlinearLeastSquares problem and use GaussNewton? That would effectively be the same algorithm as before.

7

That is surprising to me, honestly, because the NewtonRaphson method requires the jacobian matrix. Imagine the following equation:

x^2+y^2-1 = 0,

the solution curve is a unit circle. If you give only that equation to the NewtonRaphson’s method, what would you expect it to gve you? Just one point of that entire curve of solutions? And if it did that, why that point?

I do not know what NewtonRaphson() did in the previous versions, but it probably solved the nonlinear equation in a different way. The developers of this library could probably help you further.

8

Thank you for the hint. The following code runs without the error. However, it returs [NaN, NaN] as solution.

    f(u, p) = [(u[1] - p[1])^2 + (u[2] - p[2])^2]
    u0 = rand(2)
    p = [1.0, 2.0]
    prob = NonlinearLeastSquaresProblem(f, u0, p)
    sol = solve(prob, GaussNewton())

Edit:
The Levenberg-Marquardt method works for my MWE:

f(u, p) = [(u[1] - p[1])^2 + (u[2] - p[2])^2]
u0 = rand(2)
p = [1.0, 2.0]
prob = NonlinearLeastSquaresProblem(f, u0, p)
sol = solve(prob, LevenbergMarquardt(), abstol=1e-8)

@ChrisRackauckas and @RayleighLord: Thank you for your answers. I consider my problem solved.

9

Your example does not have a unique solution. I would expect that the algorithm finds one of the soultions. Which one it finds will depend on the initial guess. However, this situation can also occur with a square Jacobian.

10

@avikpal can you look into the NaNs here?

11

The condition number is blowing up due to the J^T J thing. I will fix it soon.

12

GaussNewton shouldn’t use JTJ at all, just QR?

13

Yes, the patch should be done in an hr or so