I am trying to use SciPy's scipy.optimize.least_squares to solve a system of 35 non linear equations related to multilateration.

For context :

There are 5 "reference points" which are fixed in space, but have unknown coordinates. Then there are 7 "measured points", which are also fixed in space but also have unknown coordinates. The goal is to find the coordinates of the 7 measured points in a local reference system (placed on one of the reference points). The distances between the measured points and the reference points are measured using a high precision tool (nanometer scale). However, this distance is not actually measured between the measured point and the reference point, but rather between a non-existent point between the measured point and the reference point, and the reference point. The distance between the measured point and the non-existent point is called the "deadzone", and the deadzone is the same for all distances measured, but has an unknown length.

We can summarize this problem as : trying to find the intersections of 5 spheres in 3D space 7 times (in different configurations), by only knowing the radii of the spheres, but with the added difficulty of the deadzone.

The distance between a measured point and a reference point can be written as such:

where there is an i reference point (x_i ; y_i ; z_i) and a j measured point (x_j ; y_j ; z_j). D_ij is the measured distance and Dz is the deadzone.

To find the coordinates of the reference points, the coordinates of the measured points, and the value of deadzone, I am using SciPy's least squares optimization function. The goal is to be able to feed the program approximate values for all coordinates and the deadzone (at 1e0 precision) and for the program to return optimized coordinates (1e-6 precision). It does this by trying to make the difference between the measured distance and the computed distance (distance obtained from the coordinates, which are at first the initial values).

The problem :

I now have a program that is able to run well enough to return approximate solutions, and succeeds in sanity tests, but is unable to return values that are precise enough (1e-1 precision). I tested this program by taking one of the "real" values (obtained using simulation in 3D CAD software) for a coordinate of a measured point and adding a bit of noise to it (z_7 = 28,52317743 + 1). I have tried to modify the parameters from SciPy's function, and have also tried to integrate mpmath, which is unsupported by SciPy's least squares, into the program to attempt to reach a higher level of precision. I also tried to use a sum of all residuals for optimization, rather than using the residual from a single equation to optimize only that equation. This does give better results, but they are still not precise enough.

Part of the issue is that I am still a novice in programming, and don't understand python very well. I am using AI to code, but this is not good enough.

Here is the code to compute the "residual" (the difference between the measured distance and the computed distance) :

    # Calculate the residuals
    residuals = []
    for j, (xj, yj, zj) in enumerate(measured_points):
        
        residual = (
            np.sqrt(xj**2 + yj**2 + zj**2) - (dA[j] + Dzone) +
            np.sqrt((xj - xB)**2 + yj**2 + zj**2) - (dB[j] + Dzone) +
            np.sqrt((xj - xC)**2 + (yj - yC)**2 + zj**2) - (dC[j] + Dzone) +
            np.sqrt((xj - xD)**2 + (yj - yD)**2 + (zj - zD)**2) - (dD[j] + Dzone) +
            np.sqrt((xj - xE)**2 + (yj - yE)**2 + (zj - zE)**2) - (dE[j] + Dzone)
        )**2
        residuals.append(residual)

    return np.array(residuals)

Here is the code that uses the least squares function :

    # Perform non-linear least squares regression
    result = least_squares(
        objective_function, initial_params, jac='3-point', diff_step=1e-8,
        ftol=1e-15, xtol=1e-15, gtol=1e-15, x_scale=1, method='dogbox',
        loss='cauchy', tr_solver='lsmr', max_nfev=100000000, args=(distances,), verbose=2
    ) # Parameters used here were tested to try to get better precision

    # Optimal values
    optimal_params = result.x

I know that the precision I am looking for is possible, because I have access to a program on Maple that does this computation using the same equations and using least squares regression and it is able to reach precision levels o up to 1e-12 when noise is added to a single variable.

Am I using the least squares function correctly ? What are some ways i could try to get higher precision results with python ?

I hope this describes the problem well enough, though I can provide more info if necessary. Since this is my first time asking a question on a forum, it's possible I fot something. Please do not hesitate to mention it if that is the case.

Thanks in advance for any help given.

I am trying to use SciPy's scipy.optimize.least_squares to solve a system of 35 non linear equations related to multilateration.

For context :

There are 5 "reference points" which are fixed in space, but have unknown coordinates. Then there are 7 "measured points", which are also fixed in space but also have unknown coordinates. The goal is to find the coordinates of the 7 measured points in a local reference system (placed on one of the reference points). The distances between the measured points and the reference points are measured using a high precision tool (nanometer scale). However, this distance is not actually measured between the measured point and the reference point, but rather between a non-existent point between the measured point and the reference point, and the reference point. The distance between the measured point and the non-existent point is called the "deadzone", and the deadzone is the same for all distances measured, but has an unknown length.

We can summarize this problem as : trying to find the intersections of 5 spheres in 3D space 7 times (in different configurations), by only knowing the radii of the spheres, but with the added difficulty of the deadzone.

The distance between a measured point and a reference point can be written as such:

where there is an i reference point (x_i ; y_i ; z_i) and a j measured point (x_j ; y_j ; z_j). D_ij is the measured distance and Dz is the deadzone.

To find the coordinates of the reference points, the coordinates of the measured points, and the value of deadzone, I am using SciPy's least squares optimization function. The goal is to be able to feed the program approximate values for all coordinates and the deadzone (at 1e0 precision) and for the program to return optimized coordinates (1e-6 precision). It does this by trying to make the difference between the measured distance and the computed distance (distance obtained from the coordinates, which are at first the initial values).

The problem :

I now have a program that is able to run well enough to return approximate solutions, and succeeds in sanity tests, but is unable to return values that are precise enough (1e-1 precision). I tested this program by taking one of the "real" values (obtained using simulation in 3D CAD software) for a coordinate of a measured point and adding a bit of noise to it (z_7 = 28,52317743 + 1). I have tried to modify the parameters from SciPy's function, and have also tried to integrate mpmath, which is unsupported by SciPy's least squares, into the program to attempt to reach a higher level of precision. I also tried to use a sum of all residuals for optimization, rather than using the residual from a single equation to optimize only that equation. This does give better results, but they are still not precise enough.

Part of the issue is that I am still a novice in programming, and don't understand python very well. I am using AI to code, but this is not good enough.

Here is the code to compute the "residual" (the difference between the measured distance and the computed distance) :

    # Calculate the residuals
    residuals = []
    for j, (xj, yj, zj) in enumerate(measured_points):
        
        residual = (
            np.sqrt(xj**2 + yj**2 + zj**2) - (dA[j] + Dzone) +
            np.sqrt((xj - xB)**2 + yj**2 + zj**2) - (dB[j] + Dzone) +
            np.sqrt((xj - xC)**2 + (yj - yC)**2 + zj**2) - (dC[j] + Dzone) +
            np.sqrt((xj - xD)**2 + (yj - yD)**2 + (zj - zD)**2) - (dD[j] + Dzone) +
            np.sqrt((xj - xE)**2 + (yj - yE)**2 + (zj - zE)**2) - (dE[j] + Dzone)
        )**2
        residuals.append(residual)

    return np.array(residuals)

Here is the code that uses the least squares function :

    # Perform non-linear least squares regression
    result = least_squares(
        objective_function, initial_params, jac='3-point', diff_step=1e-8,
        ftol=1e-15, xtol=1e-15, gtol=1e-15, x_scale=1, method='dogbox',
        loss='cauchy', tr_solver='lsmr', max_nfev=100000000, args=(distances,), verbose=2
    ) # Parameters used here were tested to try to get better precision

    # Optimal values
    optimal_params = result.x

I know that the precision I am looking for is possible, because I have access to a program on Maple that does this computation using the same equations and using least squares regression and it is able to reach precision levels o up to 1e-12 when noise is added to a single variable.

Am I using the least squares function correctly ? What are some ways i could try to get higher precision results with python ?

I hope this describes the problem well enough, though I can provide more info if necessary. Since this is my first time asking a question on a forum, it's possible I fot something. Please do not hesitate to mention it if that is the case.

Thanks in advance for any help given.

• python
• scipy
• precision
• scipy-optimize

• It's clear that you are new to programming. Your core problem here is that you're running before you can walk. If you can't solve y=ax+b with least_squares( ) yet, you shouldn't try to solve a set of 35 non-linear equations. This site requires a minimal-but-complete example to show your problem, but I bet you cannot minimize your problem to y=ax+b. – MSalters Commented Apr 1 at 13:27
• don't understand python very well [...] I am using AI to code - Using AI is not a very good path to learning. – Reinderien Commented Apr 1 at 23:17

1 Answer 1

You can't use higher than 64-bit floats inside least_squares - the assumption that you are using 64-bit floats or smaller is baked pretty deeply into that code.

However, your function or Jacobian can use larger than 64-bit float precision internally, and then convert the final result to a NumPy array of dtype np.float64.

Why does this help you?

One of the steps that least_squares() is doing that costs it a lot of precision is numerically differentiating your function. If you have the function f(x), then numerically differentiating it with the formula (f(x + h) - f(x)) / h can be imprecise if f(x) is larger than f'(x), which it usually is. (The differentiation method you have selected, 3-point, is not doing exactly this, but this is still a helpful intuition.)

For example, when I numerically differentiate the function f(x) = x ** 3 - 100, and compare the result to the analytical result from f'(x) = 3 * x ** 2, the result is typically only accurate to 9 digits over the range [1, 100] rather than the 15-ish that a 64-bit float is capable of.

For this reason, one idea would be to write a function which numerically differentiates the function and finds its Jacobian, using higher-precision math and a very small eps. Since this is the most imprecise step, I think this would give you the largest benefit from mpmath, and it's something you can do with the existing API of least_squares().