NLsolve.jl

Solving non-linear systems of equations in Julia.

NLsolve.jl is part of the JuliaNLSolvers family.

Non-linear systems of equations

The NLsolve package solves systems of nonlinear equations. Formally, if F is a multivalued function, then this package looks for some vector x that satisfies F(x)=0 to some accuracy.

The package is also able to solve mixed complementarity problems, which are similar to systems of nonlinear equations, except that the equality to zero is allowed to become an inequality if some boundary condition is satisfied. See further below for a formal definition and the related commands.

There is also an identical API for solving fixed points (i.e., taking as input a function F(x), and solving F(x) = x).

Note, if a single equation and not a system is to be solved and performance is not critical, consider using Roots.jl. If you want to solve a small system of equations at high performance, consider using NonlinearSolve.jl.

A super simple example

We consider the following bivariate function of two variables:

(x, y) -> [(x+3)*(y^3-7)+18, sin(y*exp(x)-1)]

In order to find a zero of this function and display it, you would write the following program:

using NLsolve

function f(x)
    [(x[1]+3)*(x[2]^3-7)+18,
    sin(x[2]*exp(x[1])-1)]
end

sol = nlsolve(f, [ 0.1, 1.2])
sol.zero

The first argument to nlsolve is the function to be solved which takes a vector as input and returns the residual as a vector. The second argument is the starting guess for algorithm. The sol.zero retrieves the solution, if converged.

A simple example

Continuing on the same system of equations, but now using an in-place function and a user-specified Jacobian for improved performance:

using NLsolve

function f!(F, x)
    F[1] = (x[1]+3)*(x[2]^3-7)+18
    F[2] = sin(x[2]*exp(x[1])-1)
end

function j!(J, x)
    J[1, 1] = x[2]^3-7
    J[1, 2] = 3*x[2]^2*(x[1]+3)
    u = exp(x[1])*cos(x[2]*exp(x[1])-1)
    J[2, 1] = x[2]*u
    J[2, 2] = u
end

nlsolve(f!, j!, [ 0.1; 1.2])

First, note that the function f! computes the residuals of the nonlinear system, and stores them in a preallocated vector passed as first argument. Similarly, the function j! computes the Jacobian of the system and stores it in a preallocated matrix passed as first argument. Residuals and Jacobian functions can take different shapes, see below.

Second, the nlsolve function returns an object of type SolverResults. In particular, the field zero of that structure contains the solution if convergence has occurred. If r is an object of type SolverResults, then converged(r) indicates if convergence has occurred.

Ways to specify the function and its Jacobian

There are various ways of specifying the residuals function and possibly its Jacobian.

With functions modifying arguments in-place

This is the most efficient method, because it minimizes the memory allocations.

In the following, it is assumed that you have defined a function f!(F::AbstractVector, x::AbstractVector) or, more generally, f!(F::AbstractArray, x::AbstractArray) computing the residual of the system at point x and putting it into the F argument.

In turn, there 3 ways of specifying how the Jacobian should be computed:

Finite differencing

If you do not have a function that compute the Jacobian, it is possible to have it computed by finite difference. In that case, the syntax is simply:

nlsolve(f!, initial_x)

Alternatively, you can construct an object of type OnceDifferentiable and pass it to nlsolve, as in:

initial_x = ...
initial_F = similar(initial_x)
df = OnceDifferentiable(f!, initial_x, initial_F)
nlsolve(df, initial_x)

Notice, we passed initial_x and initial_F to the constructor for df. This does not need to be the actual initial x and the residual vector at x, but it is used to initialize cache variables in df, so the types and dimensions of them have to be as if they were.

Automatic differentiation

Another option if you do not have a function computing the Jacobian is to use automatic differentiation, thanks to the ForwardDiff package. The syntax is simply:

nlsolve(f!, initial_x, autodiff = :forward)

Jacobian available

If, in addition to f!(F::AbstractArray, x::AbstractArray), you have a function j!(J::AbstractArray, x::AbstractArray) for computing the Jacobian of the system, then the syntax is, as in the example above:

nlsolve(f!, j!, initial_x)

Again it is also possible to specify two functions f!(F::AbstractArray, x::AbstractArray) and j!(J::AbstractArray, x::AbstractArray) that work on arbitrary arrays x.

Note, that you should not assume that the Jacobian J passed into j! is initialized to a zero matrix. You must set all the elements of the matrix in the function j!.

Alternatively, you can construct an object of type OnceDifferentiable and pass it to nlsolve, as in:

df = OnceDifferentiable(f!, j!, initial_x, initial_F)
nlsolve(df, initial_x)

Optimization of simultaneous residuals and Jacobian

If, in addition to f! and j!, you have a function fj!(F::AbstractArray, J::AbstractArray, x::AbstractArray) that computes both the residual and the Jacobian at the same time, you can use the following syntax

df = OnceDifferentiable(f!, j!, fj!, initial_x, initial_F)
nlsolve(df, initial_x)

If the function fj! uses some optimization that make it cost less than calling f! and j! successively, then this syntax can possibly improve the performance.

Providing only fj!

If a function is available for calculating residuals and the Jacobian, there is a special syntax for an, arguably, simpler approach. First, define the function as

function myfun!(F, J, x)
    # shared calculations begin
    # ...
    # shared calculation end
    if !(F == nothing)
        # mutating calculations specific to f! goes here
    end
    if !(J == nothing)
        # mutating calculations specific to j! goes
    end
end

and solve using

nlsolve(only_fj!(myfun), initial_x)

This will make enable nlsolve to efficiently calculate F(x) and J(x) together, but still be efficient when calculating either F(x) or J(x) separately.

With functions returning residuals and Jacobian as output

Here it is assumed that you have a function f(x::AbstractArray) that returns a newly-allocated vector containing the residuals. Simply pass it to nlsolve, and it will automatically detect if f is defined for one or two arguments:

nlsolve(f, initial_x)

Note, that this means that if you have a function f with a method that accepts one argument, and another method that accepts two arguments, it will assume that the two argument version is a mutating f, such as described above.

Via the autodiff keyword both finite-differencing and autodifferentiation can be used to compute the Jacobian in that case.

If, in addition to f(x::AbstractArray), there is a function j(x::AbstractArray) returning a newly-allocated matrix containing the Jacobian, we again simply pass these to nlsolve:

nlsolve(f, j, initial_x)

If, in addition to f and j, there is a function fj returning a tuple of a newly-allocated vector of residuals and a newly-allocated matrix of the Jacobian, the approach is the same:

nlsolve(f, j, fj, initial_x)

With functions taking several scalar arguments

If you have a function f(x::Float64, y::Float64, ...) that takes the point of interest as several scalars and returns a vector or a tuple containing the residuals, you can use the helper function n_ary. The complete syntax is therefore:

nlsolve(n_ary(f), initial_x)

Finite-differencing is used to compute the Jacobian.

If the Jacobian is sparse

If the Jacobian of your function is sparse, it is possible to ask the routines to manipulate sparse matrices instead of full ones, in order to increase performance on large systems. This means that we must necessarily provide an appropriate Jacobian type so the solver knows what to feed j!.

df = OnceDifferentiable(f!, j!, x0, F0, J0)
nlsolve(df, initial_x)

It is possible to give an optional third function fj! to the constructor, as for the full Jacobian case.

Note that the Jacobian matrix is not reset across function calls. As a result, you need to be careful and ensure that you don't forget to overwrite all nonzeros elements that could have been initialized by a previous function call. If in doubt, you can clear the sparse matrix at the beginning of the function. If J is the sparse Jacobian, this can be achieved with:

fill!(a, 0)
dropzeros!(a) # if you also want to remove the sparsity pattern

Fine tunings

Three algorithms are currently available. The choice between these is achieved by setting the optional method argument of nlsolve. The default algorithm is the trust region method.

Trust region method

This is the well-known solution method which relies on a quadratic approximation of the least-squares objective, considered to be valid over a compact region centered around the current iterate.

This method is selected with method = :trust_region.

This method accepts the following custom parameters:

• factor: determines the size of the initial trust region. This size is set to the product of factor and the euclidean norm of initial_x if nonzero, or else to factor itself. Default: 1.0.
• autoscale: if true, then the variables will be automatically rescaled. The scaling factors a