Nonlinear System of Equations with Bounds/Constraints on Unknowns

@bielim Glad to help! KNITRO is a mature solver that contains many people’s many years of research and efforts. It also contains several algorithms (some of the interior-point variety, including one that is quite similar to IPOPT, and some of the active-set variety). If you don’t specify an algorithm, KNITRO chooses one for you automatically at the beginning (you should see this if you inspect the output). You’ll probably want to try each algorithm and identify the one that’s best for your problem.

Using our interface to KNITRO, you can specify a starting point with

stats = knitro(model, x0=[0.5; 0.5; 0.5; 0.5; 0.5; 0.5])

or

stats = knitro(model, x0=model.meta.x0)  # use starting point specified in the model

By default, we decided to let KNITRO compute its own. The interior-point algorithms will reject your starting point if it doesn’t strictly satisfy the bound constraints.

I’m quite glad to see that TRON isn’t doing too badly either (at least with the NLS model). It’s a projected-direction method. I would have to investigate a little to determine why it has trouble solving the other model.

Feel free to open issues on our repos if you have difficulties. Also feel free to ask questions on our Slack channel: https://optimizers-workspace.slack.com/

1 Like

Try running the tuner to see if it helps

m = Model(with_optimizer(KNITRO.Optimizer, tuner = 1))

It might find an even better algorithm

This is beautiful, @saschatimme! I’ve been trying the approach you posted with different values of the parameter M – for a single parameter, homotopy continuation is slower than the other methods I tried, but the call to monodromy_solve is a one-off cost that I’m happy to pay given that I’ll have to repeat the calculation for many different values of M (and the subsequent iterations are pretty fast!).

As an aside, I’d like to express my appreciation for the amazing Nextjournal documentation – really neat and user-friendly!

When trying different values for M, I always got four real solutions. Two of them are actually the same, because the solution is symmetric in the sense that the two components of the mixture distribution the solution describes can be relabeled while leaving the mixture pdf invariant. I’m only interested in the mixture pdf, so the “double presence” of solutions is not really a problem - picking either one of the solutions is fine. (But I saw that HomotopyContinuation’s GroupActions would probably be able to filter the solutions down to the ones that are truly different, and I will give this a try.)

However, I also tried to incorporate the equality constraints on the weights (w_1+w_2= 1) into the system of equations by replacing all instances of w_2 with 1-w_1, and when I do that (now dealing with a system of 6 equations for 5 unknowns), the solver doesn’t find a real solution anymore.
I guess I could just pick the solution that satisfies the constraint on the weights most accurately, but I was wondering why directly incorporating it into the system doesn’t seem to work, and if there is a way I can include the constraint in the problem setup rather than after the solution has been computed.

Happy to hear that you find this approach helpful and that the documentation is helpful as well :slight_smile:

When you have a more equations than unknowns, then you need a very special set of equations such that you still have solutions. Geometrically, your 6 equations in 6 unkowns already only result in 18 points, and now you try to intersect these with the hyperplane 1=w_1 + w_2. So these 18 points have to be quite special such that at least one lies on this plane. Obviously, you know for theoretical reasons that this is the case here. However, I assume your measure are not exact, right?
This perturbation in the measure then would yield solutions which are a little bit away from the hyperplane. When you look at the 4 real solutions for your example above, then you notice that w_1 and w_2 don’t sum up to 1 exactly

julia> 0.31831615942330055 + 0.6949587215517059
1.0132748809750065

julia> 0.6421990255933143 + 0.38742674206054184
1.0296257676538563

If this constraint is important to be satisfied, maybe you could use the obtained solutions coming from HomotopyContinuation.jl as a starting point for a least square optimization routine.

2 Likes

many people’s many years of research and efforts

So no magic bullet, just the magic of hard work :slight_smile:

I see - I had assumed that the starting value x0 given as input to the Model (e.g., model = ADNLSModel(c, x0, 7, lvar=zeros(6))) would automatically be used as the starting point for knitro as well.

I’ll keep experimenting and will post any issues I encounter!

m = Model(with_optimizer(KNITRO.Optimizer, tuner = 1))

Apparently with_optimizer has been deprecated and replaced by Model(optimizer_with_attributes(KNITRO.Optimizer, "tuner" => 1))

I tried this – the Knitro Tuner decided that the fastest solution would be found by a solver that involves “algorithm 3”, which based on the user manual is an Active-Set algorithm. But the tuning did not lead to an improvement in the time to solution.

Yeah, you never know. Your problem is sufficiently small that there may not be much difference.

The other two very useful options for knitro to keep in mind in the future are: (1) ms_enable = 1 which turns on multi-start, and is huge time save if you have ugly, non-convex equations; and (2) honorbnds = 1 which was necessary for me when my box-bounds defined tight constraints on when my equations could be evaluated.

That makes sense, thanks for the explanation! Yes, as you point out, for theoretical reasons this overdetermined system should have a unique solution (up to flipping the labels of the two mixture components), but I understand that numerical inaccuracies can stand in the way of finding that solution.

I’ll test if using the homotopy continuation solution as an extremely educated inital guess for a subsequent optimization algorithm is too computationally expensive. As a simpler alternative, picking the solution that best satisfies the constraint, followed maybe by “re-normalizing” the weights (w_{1, \textrm{normalized}} = \frac{w_1}{w_1+w_2}, w_{2, \textrm{normalized}} = \frac{w_2}{w_1+w_2}) such that the constraint is exactly satisfied may already be good enough an approximation for my purpose.

On a different note: Is it possible to accept more than one reply as solutions? I tried to mark both @dpo’s and @saschatimme’s replies, but when marking the second one the first one gets unmarked.

Unfortunately, it isn’t possible. Just mark the one you prefer.

1 Like