# Minimization of non-linear function

Problem Summary
I am trying to estimate the (x,y) coordinates of each node in a graph, where I know the distances between connected nodes. For example

Given this graph, I can define an equation for the length of each link simply using the Pythagorean theorem.

$$L_{AB}=\sqrt{(X_B-X_A)^2+(Y_B-Y_A)^2}$$ $$L_{AD}=\sqrt{(X_D-X_A)^2+(Y_D-Y_A)^2}$$ $$L_{BC}=\sqrt{(X_C-X_B)^2+(Y_C-Y_B)^2}$$ $$L_{BD}=\sqrt{(X_D-X_B)^2+(Y_D-Y_B)^2}$$ $$L_{CD}=\sqrt{(X_D-X_C)^2+(Y_D-Y_C)^2}$$

Solution Approach

In this problem I know all of the L values, and I need to solve for all the X and Y values. I can arbitrarily pick a node, and call it the origin, so I can say A = (0,0), then I have 6 unknowns $X_{B-D}$ and $Y_{B-D}$.

I figure I can pose this question as minimization of the sum of the square of the errors, mathematically

$$f=(\sqrt{(X_B-X_A)^2+(Y_B-Y_A)^2}-L_{AB})^2\\+(\sqrt{(X_D-X_A)^2+(Y_D-Y_A)^2}-L_{AD})^2\\+(\sqrt{(X_C-X_B)^2+(Y_C-Y_B)^2}-L_{BC})^2\\+(\sqrt{(X_D-X_B)^2+(Y_D-Y_B)^2}-L_{BD})^2\\+(\sqrt{(X_D-X_C)^2+(Y_D-Y_C)^2}-L_{CD})^2)$$

Question
I'm not sure how to perform a minimization of a function in C++ (I'm not too familiar with this area, so I apologize if that is a broad question). If I were in MATLAB, I could make some initial guesses

X = [0 1 2 3];
Y = [0 2 1 3];


Then call fminunc

fminunc(@f, [X,Y])


This would return an array that I could pull apart into the X and Y coordinates.

Basically in C++ I can define some function

double f(double* A)
{
// Let A be an array of all X and Y coordinates
// This function is the error function described above
}


I would like to write (or hopefully find!) a function that does the same thing as fminunc but in C++, so I guess it would look like

double* fminunc(double(*f)(double*), A);


Where f is the pointer to my function to minimize, and A is my array of initial guesses, and it would return an array of the actual solution values.

Does such a function exist already? Otherwise, is there a name of a technique I can search for, again I'm not too familiar with this field, so something like "Jones' Method" could get me started in the right direction, I just need to know what I'm looking for.

Edit
As @TylerOlsen indicated, depending on the number of nodes and number of links, it is possible that this system may be under-constrained. Are there techniques that can provide a valid solution, even if there is no unique solution? The background to this problem is I'm working on a certain flavor of shortest-path algorithm, and as long as the coordinates respect the edge lengths, I think any valid solution should work.

• Your problem is under-constrained. You have 5 equations (one for each known length) and 6 degrees of freedom. You need to add a constraint in order to obtain a unique solution. The way I'd do it is to fix the location of a second point (eg, point B) based on its distance from A. This will remove 2 degrees of freedom and 1 equation, leaving you with 4 unknowns and 4 equations. You still need to solve them, but you will at least be able to do so now. – Tyler Olsen Jul 6 '15 at 14:30
• In the example that you've drawn, the point $C$ has two possible locations (another one below $BD$). So if you were to impose an additional constraint on $AC$, then only one of those locations would be feasible. – Kirill Jul 6 '15 at 21:50
• If @CoryKramer would use the minimization approach, the fact that the problem is under-constrained results only in having multiple local optima that have the same objective function value. If it doesn't matter for the application which minimum is found, one doesn't need to add some extra constraint. – GertVdE Jul 7 '15 at 4:58
• And instead of fixing another point, a more general approach for regularization could be to add to the objective function a measure for the size of the total graph: the distance between the origin (the fixed point) and the farthest point, i.e. $max_j \sqrt{X_j^2 + Y_j^2}$ (if needed with a weight $\theta$ in order to make it more/less important than the other constraints). – GertVdE Jul 7 '15 at 5:02

I think better approach is to use 'fzero()' function like this: then simply use (1) and (2) to find Y_B and Y_D

use same trick for other triangles

• That is also an interesting approach, I will look into that. Though as far as using fzero, again I'm working on a C++ solution so I will not have MATLAB functions available. – Cory Kramer Jul 6 '15 at 21:31
• fzero use a simple algorithm and so simple to implement, see numerical recipes chapter 9 , or this lecture http://ocw.mit.edu , also you can use open libraries e.g Madagascar – Reza Afzalan Jul 6 '15 at 21:57
• if you want source check Zero finder – Reza Afzalan Jul 6 '15 at 22:06

There are many methods to minimize a multi-dimensional function. Here is a list of some standard ones.

If you want a good description along with some sample code you can check out Numerical Recipes. In particular checkout chapter 10 which covers minimization and maximization of functions. I've implemented their simplex method to solve 2D problems like yours with good success.

Although as with most non-linear solvers, it requires an initial estimate. You'll find the same issue in your problem as you have probably already noticed that there are multiple solutions since you can flip triangles and get the same solution.

After looking through various non-linear solver libraries, I've found one that adequately solves my problem, NLopt. They have implemented a variety of optimization algorithms that are fairly well documented. The library itself has a C++ API (among other available languages)

For the problem I described, the algorithm I selected was nlopt::algorithm::LN_PRAXIS. This is described as

"PRAXIS" gradient-free local optimization via the "principal-axis method" of Richard Brent, based on a C translation of Fortran code downloaded from Netlib

This algorithm was originally designed for unconstrained optimization. In NLopt, bound constraints are "implemented" in PRAXIS by the simple expedient of returning infinity (Inf) when the constraints are violated (this is done automatically—you don't have to do this in your own function). This seems to work, more-or-less, but appears to slow convergence significantly. If you have bound constraints, you are probably better off using COBYLA or BOBYQA.

Therefore, a gradient function (Jacobian) was not required, and the algorithm applied to unconstrained systems, which applied to my problem since the solution space for the unknown variables is all reals.

To use this library, here is an example of what I had to do

Set up an nlopt::opt which is the optimizer object.

nlopt::opt optimizer{ nlopt::algorithm::LN_PRAXIS, xy.size()};


I passed it my algorithm choice, and the number of unknown variables (which were in a vector xy). Then I called set_min_objective

optimizer.set_min_objective(coord_estim::ObjectiveFunctionCallback, static_cast<void*>(&coord_est));


The ObjectiveFunctionCallback was a wrapper around my function f that I described in my post, and the requirement for the callback function was that it has the signature

double ObjectiveFunctionCallback(const std::vector<double> &x, std::vector<double> &grad, void* f_data);


I then set the stop condition, in my case I specified that the total error should be below some threshold, though there are other stop criteria.

optimizer.set_stopval(0.01);


Finally, you call optimize to perform the optimization

optimizer.optimize(xy, total_error);


The first argument of optimize will contain the optimized values, and the second argument contains the final value of your objective function corresponding to these values. optimize also provides a return code so you know what caused the optimization to terminate.

• A few comments/suggestions on your approach going forward: (1) this approach is unlikely to scale beyond a few nodes, (2) assuming you can avoid issues of nonsmoothness with the square roots (i.e., derivatives don't exist when the argument of the square root is zero), a gradient-based optimizer is likely to converge faster. IPOPT is usually the first software package recommended for gradient-based optimization of nonlinear programs in C++. – Geoff Oxberry Jul 8 '15 at 6:57
• @GeoffOxberry Thanks for the feedback! I also got a reasonable answer using LD_SLSQP which is a gradient-based, local optimization method. I can test how these various algorithms scale once I start passing in real data sets. I'll definitely consider IPOPT if performance starts being a problem. I just went with NLopt because the API was very straight-forward to use. – Cory Kramer Jul 8 '15 at 13:33