[ANN] MDBM.jl - A root-finding package for system of non-linear equations

Package Announcements
1

Multi-Dimensional Bisection Method

MDBM.jl is a package is an efficient and robust root-finding algorithm, which can be used to determine whole high-dimensional submanifolds (points, curves, surfaces…) of the roots of implicit non-linear equation systems, especially of the cases, where the number of unknowns surpasses the number of equations.

It uses far less function evaluation than the brute-force approach, making it much faster and more memory efficient, especially for complex tasks.

MDBM algorithm can handle:

  • multiple solutions
  • arbitrary number of variable (typically: 3-6)
  • arbitrary number of implicit equations
  • multiple constraints in the parameter space
  • degenerated functions

while providing the gradients of the equations for the roots by means of first order interpolation (and convergence rate).

This method can be used to create contour plots or isosurfaces in higher dimension, furthermore, it has as the advantage of being able to handle multiple functions at once.

It is a Julia implementation of a Matlab package.
Please have a look at MDBM.jl and at the examples. Feel free to suggest improvements, design choices and submit PR.

I hope you will find the algorithms useful in your own work.

Best wishes,
Daniel Bachrathy

21 Likes
2

Neat! Can you elaborate a bit about the requirements for the functions it can handle? Eg I understand that they need be continuous, but do they need to be ADable or evaluating reals to reals is sufficient?

Also, can your provide some intuition about how the solution time scales with the dimensionality of the unknowns and the number of equations?

3

Very cool! In addition to Tamas’ questions would you mind giving a bit of a run-down of how this compares to NLsolve.jl, ie. relative strengths and weaknesses?

2 Likes
4

There is no any restriction for the function, because actually the method finds the ‘sing-change’ and not the ‘root’.
So in case of f=1/x, the solution provided by the method will converge to 0 (despite to the ill conditioned interpolation).

Even, degenerated functions can be used (see the example at the end)

So it is sufficient that f evaluates reals^n to reals^m (if m<=n).

The complexity of the method is proportional to the (fractal) dimension of the object
e.g:
3 param & 2 equation → curves: 1D object, so the computation time doubles (2^1) every iteration steps
4 param & 2 equation → surface: 2D object, so the computation time quadruple every iteration steps

I forgot to include, but I gave a talk at JuliaCon2018/MDBM

The computation needs are explained at 8:09.

Degenerated example:

ax1=Axis(-3.0:3.0,"x")
ax2=Axis(-3.0:3.0,"b")

function foo(x,y)
    if norm([x,y])<1.0
        0.0
    else
        x-y
    end
end

mymdbm=MDBM_Problem(foo,[ax1,ax2])#ax3
iteration=5 #number of refinements (resolution doubling)
solve!(mymdbm,iteration,interpolationorder=1);

#evaluated points
x_eval,y_eval=getevaluatedpoints(mymdbm)

#solution points
x_sol,y_sol=getinterpolatedsolution(mymdbm)

fig = figure(1);clf()
scatter(x_eval,y_eval,s=2)
scatter(x_sol,y_sol,s=4);

degen

1 Like
5

I did not used the NLsolve.jl, however, it seems that is provides a single solution only.
The great strengths of my method that is can solve problem with more parameter than equations
e.g.:
f(x,y)=x^2+y^2-1.0
which is typical in many problem (computation of a chart, a bifurcation curve,…).
Furthermore, it can handle non-smooth functions (see the previous comment).
This type of problem, where the existence of multiple solution is important (especially as a function of one or two parameter) the MDBM is excellent.

Compared it to a continuation method, the MDBM is automatic (e.g.: no need to start a new branch), and it can find closed curves (isola) which is hard to detect by continuation.

The drawbacks of MDBM

  • that is has only linear convergence (but it is not a big problem in chart plotting)
  • it can also miss small closed curves (isola) if the initial resolution is not fine enough
  • handling more the 6 parameter or computing higher order objects (surfaces, volumes…) can be very memory & time consuming
2 Likes
6

p.s.: we have found a small bug, which is fixed in the new 0.1.3 release!

7

Thank you for creating this module. It is much appreciated.

8

Very nice! How does this compare to IntervalRootFinding.jl?

1 Like
9

By the way, I think it would be useful to change the package name to be more explicit and discoverable.

3 Likes
10

Very cool! Look forward to trying it out on some problems I have!

11

I think that the IntervalRootFinding.jl is definitely better for simpler analytical expressions than the MDBM.jl, because it can provide guaranteed and unique solutions, while in the MDBM there is slight a chance to lose some solutions, however the usage of the IntervalRootFinding.jl is limited.

The MDBM is initially designed for stability chart computations which is based on the eigenvalues of a large matrix, for which the IntervalRootFinding would fail (see the first example below).
So the MDBM can be used for broader problems, e.g.: eigenvalue based (1st example) and iteration based (2nd example) functions.
Furthermore, it was used to ‘black-box’ functions where the output of a Finite Element Software (Abaqus) is used. The MDBM method was successfully used in automatically evaluated measurements to find the stability boundary of a controlled mechanical system, where noise was significant (see the third example with nois.

1# Egienvalue problem

using LinearAlgebra, Arpack
function fooeig(x)
       mymx=Matrix{Float64}(I,5,5);
       mymx[2,1]=x;
       mymx[1,2]=x+5;
       abs(eigs(mymx)[1][1])-3
end
eig_problem=MDBM_Problem(fooeig,[Axis(-10:10)])
solve!(eig_problem,10)
getinterpolatedsolution(eig_problem)

Roots:

1-element Array{Array{Float64,1},1}:
 [-5.70156, 0.701562]

#2. Iteration based problem:

function mandelbrot(x,y)    
    c=x+y*im
    z=Complex(zero(c))
    k=0
    maxiteration=1000
    while (k<maxiteration && abs(z)<4)
            z=z^2+c
            k=k+1
        end
    return abs(z)-2
end

Mandelbrotmdbm=MDBM_Problem(mandelbrot,[-5:2,-2:2])
solve!(Mandelbrotmdbm,9)
a_sol,b_sol=getinterpolatedsolution(Mandelbrotmdbm)
fig = figure(3);clf()
plot(a_sol,b_sol,linestyle="", marker=".", markersize=1)

Mandelbrot

#3. Function with noise

function foo(x,y)
    x^2.0+y^2.0-2.0^2.0+rand()
end

mymdbm=MDBM_Problem(foo,[-3.0:3.0,-3.0:3.0])
solve!(mymdbm,5)

x_sol,y_sol=getinterpolatedsolution(mymdbm)
fig = figure(1);clf()
scatter(x_sol,y_sol,s=4);

noisy_circle

1 Like
12

Nice examples, thanks!

13

Dr. Bachrathy,

From the peanut gallery with pedantic thoughts… I comment on using acronym for the package: it could have been named MultiDimensionalBisection.jl. There’s no harm in a “using MultiDimensionalBisection” at the top of a file, and actually, it’s quite helpful/informative. On the other hand, MDBM.jl made me think it’s some sort of database or … well, acronyms are not as informative as longer, spelled-out package names. This makes it harder to search for on google as well. I don’t know if it’d even be possible to change the name, etc. Regardless, I love the examples: as someone who is not particularly numerically inclined they were quite illustrative.

Congragulations,
Clark

5 Likes