Minimization of non-linear function

Question

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.

Answer 1

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

Answer 2

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.

Answer 3

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.