[Python] A more efficient and robust algorithm for minimizing an objective with many local...

Background


Background I've developed a Python script that calculates the time taken for a ray of light to propagate between two points in a complex medium. Now, I want to extend this script to solve the Time Difference of Arrival (TDoA) problem in that medium. The problem statement is as follows:

Given the known locations of N receivers (e.g., cameras, sensors, antennae) and the time of arrival of a light signal at each receiver, I aim to determine the coordinates of the signal emission. Since the time of signal emission is unknown, I only know the pairwise arrival time differences (e.g., the difference in signal arrival time between receiver A and receiver B).

What I've Done So Far


I've implemented a script to tackle this problem, as shown below:

from dataclasses import dataclass
import numpy as np
from scipy.optimize import differential_evolution
from ray_tracer import RayTracer
from random import randrange

@dataclass
class Optimizer:
ice_model: np.ndarray # Parameters used by my RayTracer class to model complex medium
antenna_coordinates: np.ndarray

def __post_init__(self):
self.ray_tracer = RayTracer(self.ice_model)

def _objective(self, coords: np.ndarray, measured_times: np.ndarray) -> float:
"""Objective function to minimize the error between measured and predicted transit times."""
x, y, z = coords
predicted_times = np.array([
self.ray_tracer.transit_time([x, y, z], [nx, ny, nz])
for nx, ny, nz in self.antenna_coordinates
])

# Compute time difference matrices
measured_time_diff = np.abs(measured_times[:, None] - measured_times)
predicted_time_diff = np.abs(predicted_times[:, None] - predicted_times)

# Calculate error matrix and fill the diagonal with zeros
error_matrix = np.abs(measured_time_diff - predicted_time_diff)
np.fill_diagonal(error_matrix, 0)

# Return the maximum error
return np.nanmax(error_matrix)

def solve(self, measured_times: np.ndarray) -> np.ndarray:
"""Solve the optimization problem to find coordinates minimizing the objective function."""
bounds = [(-1000, 1000), (-1000, 1000), (-1000, 0)]

# Use differential evolution to minimize the objective function
result = differential_evolution(
self._objective,
bounds=bounds,
args=(measured_times,),
strategy='best1bin', # Strategy for optimization
maxiter=5000, # Maximum number of iterations
tol=1e-7, # Tolerance for convergence
popsize=30, # Population size
mutation=(0.5, 1), # Mutation factor
recombination=0.7, # Recombination constant
disp=True # Display convergence messages
)

if result.success:
return result.x # Return the optimized coordinates
else:
raise ValueError("Optimization failed: " + result.message)

# Example usage
if __name__ == "__main__":
# Define the ice model
ice_model = np.array([1.78, 0.454, 0.0132])

# Initialize RayTracer with the ice model
ray_tracer = RayTracer(ice_model)

# Generate random antenna coordinates
antenna_coordinates = np.array([[randrange(-100, 0) for _ in range(3)] for _ in range(10)])

# Generate random true coordinates and measured transit times
x, y, z = [randrange(-1000, 0) for _ in range(3)]
measured_times = np.array([
np.random.normal(ray_tracer.transit_time([x, y, z], [nx, ny, nz]) + 10, 1e-12)
for nx, ny, nz in antenna_coordinates
])

# Initialize the optimizer
optimizer = Optimizer(ice_model, antenna_coordinates)

# Attempt to find optimal coordinates
try:
optimal_coordinates = optimizer.solve(measured_times)
print("Optimal Coordinates (x, y, z):", optimal_coordinates)
print("True Coordinates (x, y, z):", x, y, z)
except ValueError as e:
print(e)



Example Output:

Optimal Coordinates (x, y, z): [-384.46240934 730.54140202 -612.43716286]
True Coordinates (x, y, z): -704 -408 -611

1. 
   The _objective function computes the signal arrival time at each receiver location for a given emission coordinate guess using the RayTracer class. I construct matrices containing the pairwise arrival time differences and return the largest absolute error. The solve method uses SciPy's differential_evolution optimizer, which I selected because it was the only one that approached solving this problem.
   
   
2. 
   The main block generates example data by choosing random receiver and source coordinates, using RayTracer to compute arrival times, and passing these times (with Gaussian noise, and an offset of 10 to prove that the signal emission time is irrelevant to the algorithm) to the solver.
   

Problem & Attempted Solutions


I require a solution that is efficient with respect to execution time, as this will (ideally) ultimately find near-real-time application. Regrettably, though I've tried applying nearly all of SciPy's optimizers (carefully considering each of their strengths, weaknesses and parameters), as well as some of my own construction, I've been unable to accurately reconstruct the source location. This does not improve with additional receivers, nor with less Gaussian noise.

The objective function is highly nonlinear, and I've tried applying a variety of smoothing techniques (Gaussian smoothing, Huber loss, Savitzky–Golay filters) to no avail. No other reasonable objective function does not suffer from this issue. Consider the following plot of the objective in z produced using example data (and with x and y fixed at their true values):



The code posted above, making use of a differential evolution solver, somewhat reliably reconstructs the z-coordinate (getting to within a unit or so in roughly 2/3 of cases, and getting nowhere near the result otherwise) but reliably fails to accurately reconstruct the x and y coordinates. It occurs to me that it would perhaps be useful to write some algorithm that first solves for the z-coordinate (reducing the problem to one in two dimensions afterwards), however I have failed to write one (I've tried e.g. sampling a handful of points from a bunch of slices of z, and taking the slice with the least error by various measures, however this was slow and unreliable).

Question: Request for Alternatives


What alternative algorithms could I employ that are faster, more accurate, and more reliable? I'd prefer an algorithm with minimal tunable parameters, as I want it to be easily applicable to various values of ice_model, bounds, etc.

Continue reading...