How to impove multistart/lsqnonlin solution selection when difference between Fvals is too small to discern from rounding noise.

I am solving a 2.5D hyperbolic acoustic TDOA localisation problem in MATLAB using the lsqnonlin optimisation solver. The objective is to minimise the variable 'residuals'. As per the documentation for lsqnonlin, I am passing 'residuals' to the solver as a signed, three element vector, rather than using the L2 norm of 'residuals' or some other scalar value. The cost function that computes 'residuals' is included below.

I am solving using "multistart", which returns multiple candidate solutions. In approximately 65% of my test cases, the system returns a reasonably good, high precision location estimate. In certain edge cases however, the difference in the final function value 'Fval' between candidate solutions is on the order of 1e-30, which means that the meaningful differences between solutions' final residual values are lost in rounding errors.

In these edge cases, the multistart solver picks the wrong solution, resulting in high localisation error. In almost all of these cases, one of the other candidate solutions has in fact found the correct location, but not seleted it, because multistart doesn't consider this solution "optimal" as defined by either Fval, first order optimality, or some combination of these metrics. I have attached a figure to illustrate this point.

In order words, the solver is converging to a valid set of local minima, one of which is almost always the global minimum, but the variation in "Fval" for these local minima is so small that we are unable to tell which is the global minimum, due to the noise in Fval caused by numeric precision errors.

I have tried calculating the residuals using VPA arrays at 35 significant digits, but lsqnonlin() only supports double precision inputs. Converting the VPA residual to double obviously results in rounding to 16 significant digits, and this is not sufficient to chose the correct solution by Fval. I have tried all kinds of residual scalings and alternate cost fuctions (huber, log1p, MLE etc.), and endless hyperparameter tinkering, but everything results in the same solution selection problem, in some non-trivial proportion of test cases.

I would really appreciate it if somebody could give me some advice for solving this.

Many thanks,

Ben

function final_residuals = calcResiduals(x, problemData)

% Insert z coord into candidate location for propagation modelling
x = [x, problemData.zVal];

% Set VPA array to use 35 significant digits
digits(35)

% Calculate simulated TDOAs - use variable precision arrays
modelTdoas = vpa(zeros(size(problemData.measuredTdoas)));
for i = 1:length(problemData.rcvPairs)
    prop1 = propagationSim(x, problemData.rcvPos(problemData.rcvPairs(i,1),:), ...
        problemData.T, problemData.S, problemData.Fs, problemData.lat, problemData.lon);
    prop2 = propagationSim(x, problemData.rcvPos(problemData.rcvPairs(i,2),:), ...
        problemData.T, problemData.S, problemData.Fs, problemData.lat, problemData.lon);
    modelTdoas(i) = prop1.propDelay - prop2.propDelay;
end

% Set timing precision (in seconds)
samplingPrecision = vpa(1/problemData.Fs);
gpsPrecision = vpa(1e-8);
timingPrecision = samplingPrecision - gpsPrecision;

% Normalize TDOAs by timing precision 
normalizedModelTdoas = modelTdoas / timingPrecision;
normalizedMeasuredTdoas = vpa(problemData.measuredTdoas) / timingPrecision;

% Calculate residuals
normalized_residuals = normalizedModelTdoas - normalizedMeasuredTdoas;

% Return the final residual as a double
final_residuals = double(normalized_residuals);

end

Answer 1

[Matt J]

I have tried calculating the residuals using VPA arrays at 35 significant digits

Do you really need to keep all 35? I imagine you wouldn't when the iterations approach convergence. Perhaps discard the most significant 20 digits and convert the final 15 back to double.

[Ben]

Hi Matt, thanks for the suggestion.

I've written a function that takes the solutions returned by Multstart() and runs a post-hoc analysis, re-calculating the Fvals as Variable Precision Arrays (digits = 50), then returns the index of the solution with the smallest Fval using these new high precision values.

Your solution is more elegant, and I have tried it, but as it turns out, neither seems to solve the problem. It seems that the correct multistart solution (the one where 'X' is closest to the true source location) is not always the one with the smallest Fval, even when machine precision errors are removed from the equation. First order optimality also does not appear to be a reliable indicator of which solution to choose.

I am becoming concerned that this problem is intractible due to the limited number of receivers, the geometry of the reeciver array, and the large range over which the source signals are able to be received.

If anybody has any advice, I would be very appreciative.

[Matt J]

Well, it does sound like you have a bad cost function formulation, although I don't know the engineering application well enough to recommend how to fix it.

But a fundamental assumption of an inverse problem is that the cost function global minimum is unique and distinct, both in its Fval and in the X that attains it. You on the other hand have a cost function which barely distinguishes the desired solution X from spurious solutions at all...

[Ben]

If the cost function wasn't so simple, I'd agree, but I can't really see where it could have gone wrong.

The application is hyperbolic localisation by time difference of arrival (TDOA). There are three co-planar, time-synchronized receivers, and one signal emitter with unknown (x,y) coordinates, but whose z coordinate is known and fixed. The signal emission time is not known. I use the time that the signal is received at each receiver to calculate the time differences of arrival between each pair of receivers. In the optimisation, this is 'measuredTdoas'.

Theoretically, the location of the receiver is the one that best minimises the error between the observed 'measuredTdoas', and 'modelTdoas', which are the TDOAs that the propagation model predicts we would see at the receivers for a source at our candidate location, x.

So the cost functions I have been using are all some normalized or scaled variation of the following:

where:

• deltaT_ij is the observed time difference of arrival for receivers i & j
• deltaT hat _ij is the predicted time difference of arrival for receivers i & j at candidate point, x.

Note that the spatial coordinates and time differences are scaled to similar orders of magnitude when passed to lsqnonlin. Whenever these values need to be physically meaningful, they are "un-scaled", eg. when 'modelTdoas' is computed for candidate point 'x', the coordinates in 'x' are returned to their original physical scale to pass to the propagation model.

I have experimented with all sorts of residual scalings and normalizations for this cost function to try to improve the conditioning of the cost surface shown below. The slopes do change, but everything I have tried still results in a ravine that is near-flat in the radial direction from origin to source. Note that this means I can correctly estimate the bearing to the source with high precision in 100% of test cases, but in certain edge cases, the geometry of the source and receivers results in the ravine being too flat along its length for Fval to be useful in discriminating between possible near-optimal solutions, thereby making range to the source along this bearing ambiguous.

Note that in my current tests, there are no actual signals involved. I use a propagation model to predict delay times for a source at a given location. I test at various randomly selected locations, and compare localisation estimates with ground truth. Also note that I am attempting to solve for location using a single measurement of TDOA, not by calculating a trajectory over time using multiple measurements, and using feasible trajectories to resolve ambiguities. I have considered doing this, but would like to avoid if possible.The propagation model has been validated and tested, and I am certain that it is not the source of errors.

It is possible that I may just be reaching the limits of the source-receiver geometry under this method. There is a moderate correlation between the true source-receiver range, and L2 norm of location estimate error. This makes sense, because in the geometry that underpins this localisation method, the source lies on the intersection point between three hyperbola branches (the hyperbola geometries are defined by the TDOAs, converted to distance differenes).

For source locations whose range to the receivers is extremely large relative to the spacing of the receiver array, I'm pretty sure the hyperbola branches approach co-linear, and so intuitively, it seems that there may be no single intersection point, but rather a region of overlap where location is ambiguous.

This is a little outside my domain though, so I'm not sure if my intuitions are correct...