forum.alglib.net

Hi,
Apologies for the lengthy context before asking my question.

I have a simple/small subproblem that I will solve thousands of times within a framework to solve a larger problem; the subproblem maximizes the distance (for this purpose, Euclidean only) between two convex neighborhoods. For simplicity's sake, let them be two ellipsoids with explicitly fairly different radii with no additional linear constraints imposed.

I have multiple means of solving this problem to global optimality including, now (thank you!!!!), alglib's options (filter SQP, AUL, etc) that I just implemented/tested; I also have tested SNOPT via dll, SNOPT, CONOPT, MINOS and COUENNE via GAMS, and a geometric parameter search method I wrote. I say 'globally optimal' on purpose because I have a geometric partitioning scheme that mimics what a branch/cut implementation would do here when used in conjunction with a local optimizer. I provide more context in this recent unanswered post: https://math.stackexchange.com/questions/4917464/analytical-solution-to-a-toy-sized-but-interesting-optimization-problem. Note that I've correctly modified the object function and gradient to solve as a minimization within ALGLIB and match the expected solution.

I solve the problem in what I call unit-canonical space (the ellipsoids are transformed into axis-aligned unit spheres centered at the Cartesian origin) because this provides numerical stability over many orders of magnitude difference in radii sizes (so my ALGLIB scaling is simply a unit-vector). The KKT conditions are quite 'simple' in that I have only 6 variables (the xyz-position for each neighborhood) and 2 quadratic constraints, so I have 6 equations to evaluate: df/dx -lambda1*dg1/dx -lambda2*dg2/dx. A further simplification is that since the dg1 vector is non-zero for only ellipsoid 1 and the dg2 vector is non-zero for only ellipsoid 2, I end up with 3 equations for df-lambda1*dg1 and 3 for df-lambda2*dg2 (because the other lambda is irrelevant in each case).

Here are two ALGLIB questions that are probably simple to you but not to me:
1. I have the analytical gradient (it's linear) and the analytical Hessian (it's a constant matrix) for f and g but there seem to be no algorithms making use of the exact Hessian - e.g. filter SQP algorithms like SNOPT talk about an approximate Hessian and a modified Lagrangian to deal with large and/or highly non-linear problems; I see nothing in minnlcsetalgosqp or the other algorithms where I can use that. Is there a means of using this within ALGLIB's optimization schemes? Are there algorithms for simple problems like mine that utilize the analytical Hessian?
2. ALGLIB and SNOPT's optimal (termination status 1) objective function value matches mine but the KKT conditions that I calculate are infeasible at the 1e-3 to 1e-5 level whereas in my (slow) geometric algorithm, I can sometimes get to 1e-11 to 1e-13 even for highly ellipsoidal cases. Even when I hack the ALGLIB sqp code to force things like egap and edual to be < 1e-14 instead of relative, the solution still is KKt feasible only at ~1e-4.
Where can I extract the options for KKT multipliers that ALGLIB computes so I can see what ALGLIB says df/dx-lambda1*dg1/dx-lambda2*dg2/dx are? All I can find are the box constraint and linear constraint KKT multipliers which are not part of my current problem? What can I do in ALGLIB to find a solution feasible to the KKT conditions at the ~e-10 level?

Thank you,
-William

If it helps, I'm willing to post the solution found by alglib fsqp, the KKT conditions I derived, the multipliers and evaluations produced from the KKT point as well as my own solution. I'm convinced there is a better way to solve this very small problem and I'd like to get the KKT feasibility down to ~e-12 in general.

Hi!

Would you like to run ALGLIB with logging enabled, i.e. call alglib::trace_file("sqp","trace.log") prior to running the optimizer? I would like to take a closer look at the optimization process. The fact that the objective value shows good agreement, but KKT error is not so good, suggests that there is some kind of degeneracy here.

Sergey.

Thank you for reading my post and or the suggestion. I've attached the trace file as well as a more detailed one.
FWIW, to show how similar the solution values are, here is the alglib solution. Please excuse the many decimal places. I generally ignore Visual Studio C++ numerical places beyond 1e-13 but I value extreme precision in my problem.
Alglib filter sqp:
Unit-Canonical Ellipsoid 1: (-0.329176838420615, 0.94401746027854, 0.0217633576606146)
Unit-Canonical Ellipsoid 2: (0.930895021061298, 0.329232641233283, -0.158241359036905)
In general space, this corresponds to
General Ellipsoid 1: (-10.5967582185939, 206.347913283321, 6.77869831250489)
General Ellipsoid 2: (18.8893523034911, -49.9697297257856, 570.805604711329)
The Euclidean norm is: 620.237467405261
In code, of the 6 KKT condition equations, the worst feasibility is 9.7502106655156240e-06 and in my Excel "implementation" of my geometric convergence algorithm, the KKT feasibility evaluates as 0.0001294260728173 or roughly 1e-4.

In contrast, the point I find in my Excel "algorithm" differs from the ALGLIB solution by at most e-9 in unit-canonical space and e-8 in general space!
My solution:
Unit-Canonical Ellipsoid 1: (-0.32917684136241, 0.944017459191415, 0.0217633603207913)
Unit-Canonical Ellipsoid 2: (0.930895019622935, 0.329232643072075, -0.158241363672689)
In general space, this corresponds to
General Ellipsoid 1: (-10.5967582026675, 206.347913213863, 6.77869828010731)
General Ellipsoid 2: (18.8893522318179, -49.9697297427802, 570.805604707352)
The Euclidean norm is: 620.237467405261
In code, of the 6 KKT condition equations, the worst feasibility is 9.0949470177292824e-12 and in my Excel "implementation" of my geometric convergence algorithm, the KKT feasibility evaluates as 0.00000000003547029 or roughly 3.5e-11.

Does this give you any hints? I realize it may seem like polishing an already very nice apple, but I'd love to squeeze out the extra accuracy from alglib if at all possible...

I forgot to give you an idea of the original problem.
Ellipsoid 1: radii: (16.765, 55.8273, 1.3), origin: (-20.0, 155.9, 19.99), rotations: (1.4, 236.83383, 12.999) degrees
Ellipsoid 2: radii: (8.7, 10.1, 15.1617), origin: (20.09, -45.6, 562.94), rotations: (1.23456, 98.764523, 186.34) degrees

So a 'large' difference in radii which I like because it shows robustness of the approach. More spherical ellipsoids solve to e-13 KKT feasibility by any method I mentioned above.

I just tested in GAMS (CONOPT and SNOPT) a "post-processing" NLP that makes the KKT multiplier calculations actual variables and constraints and minimizes the infeasibility as the objective function while requiring my previous objective function value be maintained...A poor person's dual formulation in a way. I have to check for divide by zero bu

It seems to only take a couple internal iterations (I guess that's the QP solver) to achieve what I want. However, it seems I should be able to get this from a primal solution. Recognizing ALGLIB advertises the "no tuning needed" for its excellent implementations, any thoughts on internal settings that might help here? Or parts of the code I could look into different from what I reference above?

Thank you!!

Hi! Actually, I thought for several days about using arbitrary precision to further polish solution obtained in a double precision. ALGLIB has legacy version (2.x branch) supporting MPFR, probably if you use run a 128-bit Newton iteration straight from the solution returned by ALGLIB it will converge in several iterations even without merit functions, line search and other kinds of safeguards that make optimizer so difficult to implement. That legacy version does not support advanced optimization algos, but you can use its linear algebra to implement Newton method yourself.

Your approach is also a good option, though.

Thank you! Yours is also an excellent idea. I will try to implement that if it doesn't take me too long to get it running.
In the guts of the current alglib release, when I comment out the termination conditions that specify if the change in objective function is too small and if the trust region is smaller than the sqrt of machine epsilon (optimization.cpp ~line 63417 and line 63438), the solution will be perturbed until the trust region is smaller than the actual epsilon I input (rather than the sqrt of epsilon). Of course, since it's not optimizing anything really at that point, it's "random"....but perhaps you could add a "solution polishing" step under these conditions that used extra numerical precision, per your idea, or explicitly minimizes the infeasibility per an extended formulation.
I'll let you know how your idea goes once I get to it - but may be a few weeks.

Thank you for thinking about it!
-William

Just an update - I think your Newton-Raphson suggestion is superior. Although my idea is certainly robust, it's far more complicated whereas this is, so far, 1-2 iterations of Newton-Raphson to get to 1e-13 roughly...and that's simply prototyping it in Excel.
I was unable to get the ALGLIB 2.6.0 mpr fully imported (I create dlls out of all 3rd party software I use and export methods I want but I can't figure out some of the unresolved externals I'm getting)...but given the success using regular 64-bit Excel accuracy, I expect to be just fine without it. Hopefully, if I ever have to use it, I'll be able to just pull in some of the specific code that had no issues linking.
So for now, I think this should be good enough.
Now to figure out how to calculate KKT multipliers for additional linear constraints.
Many thanks!!