Newton-Raphson implementation, SymPy linsolve producing empty set?

Question

I wrote am implementation of the Newton-Raphson method (a multivariate generalization of Newton's method) to solve the system:

sqrt( (x - x_i) ** 2 + (y - y_i) ** 2 + (z - z_i) ** 2 ) - 299792458 * (t_i - t),
sqrt( (x - x_j) ** 2 + (y - y_j) ** 2 + (z - z_j) ** 2 ) - 299792458 * (t_j - t),
sqrt( (x - x_k) ** 2 + (y - y_k) ** 2 + (z - z_k) ** 2 ) - 299792458 * (t_k - t),
sqrt( (x - x_m) ** 2 + (y - y_m) ** 2 + (z - z_m) ** 2 ) - 299792458 * (t_m - t)

Where the unknowns are x, y, z and t. This system describes the solution to a localization problem, wherein the location of three observers, the velocity of some signal (in this case, c = 299792458), and the time at which the signal reached each observer is known, and we want to determine the coordinates x, y, z of the signal source (and, in doing so, we also find the time of emission, t).

I believe to understand the algorithm well. In this smaller example for solving the system:

10 * x + 3 * y * y - 3
x * x - exp(y) - 2

I attained excellent convergence, even with not-so-good initial guesses.

from dataclasses import dataclass
from sympy import *

x, y, z = symbols('x y z')


@dataclass
class Solve:

    @staticmethod
    def newton_raphson():

        F = Matrix([10 * x + 3 * y * y - 3, x * x - exp(y) - 2])
        v = Matrix([x, y])
        J = F.jacobian(v)

        xx, yy = 100000, 100000
        for i in range(100000):
            A = N(J.subs({"x": xx, "y": yy}))
            b = N(F.subs({"x": xx, "y": yy}))
            b *= -1

            update = linsolve((A, b), [x, y])

            (dx, dy) = tuple(*update)

            xx += dx
            yy += dy

            print(xx, yy)




mySolver = Solve
mySolver.newton_raphson()

My Code:

from sympy import *

x, y, z, t = symbols('x y z t')
x_i, y_i, z_i, t_i = symbols('x_i y_i z_i t_i')
x_j, y_j, z_j, t_j = symbols('x_j y_j z_j t_j')
x_k, y_k, z_k, t_k = symbols('x_k y_k z_k t_k')
x_m, y_m, z_m, t_m = symbols('x_m y_m z_m t_m')


F = Matrix([
            sqrt( (x - x_i) ** 2 + (y - y_i) ** 2 + (z - z_i) ** 2 ) - 299792458 * (t_i - t),
            sqrt( (x - x_j) ** 2 + (y - y_j) ** 2 + (z - z_j) ** 2 ) - 299792458 * (t_j - t),
            sqrt( (x - x_k) ** 2 + (y - y_k) ** 2 + (z - z_k) ** 2 ) - 299792458 * (t_k - t),
            sqrt( (x - x_m) ** 2 + (y - y_m) ** 2 + (z - z_m) ** 2 ) - 299792458 * (t_m - t)
            ]).subs({'x_i': 1, 'y_i': 3, 'z_i': 5, 't_i': 3.678593478725817e-07,
                     'x_j': 2, 'y_j': 4, 'z_j': 6, 't_j': 3.6211297837266057e-07,
                     'x_k': 3, 'y_k': 6, 'z_k': 9, 't_k': 3.5014698296079636e-07,
                     'x_m': 5, 'y_m': 10, 'z_m': 15, 't_m': 3.26330717010967e-07})
J = F.jacobian([x, y, z, t])

#print(J, F)

def compute():
    xx, yy, zz, tt = 0, 0, 0, 0
    max_iter = 1000
    for i in range(max_iter):
        A = N(J.subs({'x': xx, 'y': yy, 'z': zz, 't': tt}))
        b = N(F.subs({'x': xx, 'y': yy, 'z': zz, 't': tt})) * -1

        update = linsolve((A, b), [x, y, z, t])

        print(update)

        (dx, dy, dz, dt) = tuple(*update)

        xx += dx
        yy += dy
        zz += dz
        tt += dt

    print('(x, y, z):', xx, yy, zz)

compute()

I "simulated" a source with some random coordinates, and chose the coordinates of the observers to be:

x_i = 1
y_i = 3
z_i = 5

x_j = 2
y_j = 4
z_j = 6

x_k = 3
y_k = 6
z_k = 9

x_m = 5
y_m = 10
z_m = 15

To produce times of arrival of:

3.678593478725817e-07 3.6211297837266057e-07 3.5014698296079636e-07 3.26330717010967e-07

My Problem:

The code above is throwing the error:

    (dx, dy, dz, dt) = tuple(*update)
ValueError: not enough values to unpack (expected 4, got 0)

It seems that linsolve is producing an EmptySet, though I'm not sure why? Note that this shouldn't be an issue with the solvability of this system. It is my understanding that Newton-Raphson is the "standard" technique used in solving this problem.

Edit:

It seems that A is:

Matrix([[-0.169030850945703, -0.507092552837110, -0.845154254728517, 299792458.000000], 
        [-0.267261241912424, -0.534522483824849, -0.801783725737273, 299792458.000000], 
        [-0.267261241912424, -0.534522483824849, -0.801783725737273, 299792458.000000], 
        [-0.267261241912424, -0.534522483824849, -0.801783725737273, 299792458.000000]])

And b is:

Matrix([[104.365378313899], 
        [101.075425086493], 
        [93.7464525227794], 
        [79.1232008397505]])

So A is singular... not sure what I've done wrong?

Update:

I've made some progress. I fed the algorithm these times:

4.456421305318992e-07 4.398748057660546e-07 4.2865257523627983e-07 4.064283790801526e-07

Which correspond to the coordinates:

83 81 76

i.e., that's what I expect my NR code to produce. By changing my initial values to 100, 100, 100 (which are fairly good starting points, I think), I was able to avoid the singular matrix issue most of the time. Now, however, I'm getting nowhere near the sort of values I'm expecting. Again, I want 83 81 76, but I'm getting (where each line is a step in Newton's method):

(x, y, z): 17456269941258.7 -34912544447499.9 17456272342079.8
(x, y, z): 227660996986076. -376466056005027. 222965719150785.
(x, y, z): -1.21023496190953e+15 2.13688425095421e+15 -726879053110979.
(x, y, z): 4.90735107057989e+17 -1.41978303824902e+17 -1.26827668619485e+17
(x, y, z): -321561212.117731*t + 0.289317600845958*y - 4.1103972009836e+17 y - 1.41978303824902e+17 -1.26827668619485e+17
(x, y, z): -321561212.117731*t + 0.289317600845958*y - 4.1103972009836e+17 + 1.0*(8.01168826792966e-19*t**2*y*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**0.5 - 8.01168826792966e-19*t**2*y*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**0.5 - 0.113748591105452*t**2*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**0.5 + 0.113748591105452*t**2*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**0.5 - 7.20833962883084e-28*t*y**2*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**0.5 + 7.20833962883084e-28*t*y**2*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**0.5 + 1.64883758053215e-9*t*y*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**0.5 + 6.55493667660098e-10*t*y*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**1.0 - 1.64883758053215e-9*t*y*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**0.5 - 6.55493667660098e-10*t*y*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**1.0 + 1.33659310664116e-26*t*z - 219568708.172346*t*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**0.5 - 93065879.1023445*t*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**1.0 + 219568708.172346*t*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**0.5 + 93065879.1023445*t*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**1.0 - 3.77928683421434e-19*y**2*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**0.5 - 5.89765954818281e-19*y**2*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**1.0 + 3.77928683421434e-19*y**2*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**0.5 + 5.89765954818281e-19*y**2*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**1.0 + 0.644246732305085*y*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**0.5 + 1.00536107970607*y*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**1.0 - 0.644246732305085*y*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**0.5 - 1.00536107970607*y*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**1.0 + 1.4351559202707e-17*z - 8.38508328353584e+16*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**0.5 - 1.30851053806647e+17*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**1.0 + 8.38508328353584e+16*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**0.5 + 1.30851053806647e+17*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**1.0)/(-3.38882283462319e-35*t**2 + 3.04901852355942e-44*t*y - 8.40342129428199e-26*t + 3.27385871089648e-35*y - 5.11605818702297e-17) y - 1.41978303824902e+17 + (-0.393024115662512*t**2*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**0.5 + 0.393024115662512*t**2*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**0.5 + 3.53614770479379e-10*t*y*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**0.5 - 3.53614770479379e-10*t*y*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**0.5 - 758653768.964884*t*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**0.5 - 321561212.117731*t*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**1.0 + 758653768.964884*t*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**0.5 + 321561212.117731*t*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**1.0 + 0.185397985565395*y*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**0.5 + 0.289317600845958*y*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**1.0 - 0.185397985565395*y*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**0.5 - 0.289317600845958*y*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**1.0 - 2.89721385578571e+17*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**0.5 - 4.52116542333159e+17*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**1.0 + 2.89721385578571e+17*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**0.5 + 4.52116542333159e+17*(4.29671948342487e-19*t**2 - 7.73175697569789e-28*t*y + 1.2082414688189e-9*t + 4.50319375723489e-36*y**2 - 2.26703331823469e-18*y + 1)**1.0)/(5.16987882845642e-26*t + 5.55111512312578e-17) z - 1.26827668619485e+17

Which is... not good.

user_9450991 · Answer

If x, y and z are all equal then the matrix is singular:

In [26]: F                                                                                                                        
Out[26]: 
⎡                  ________________________________                    ⎤
⎢                 ╱        2          2          2                     ⎥
⎢ 299792458⋅t + ╲╱  (x - 1)  + (y - 3)  + (z - 5)   - 110.281458096998 ⎥
⎢                                                                      ⎥
⎢                  ________________________________                    ⎥
⎢                 ╱        2          2          2                     ⎥
⎢ 299792458⋅t + ╲╱  (x - 2)  + (y - 4)  + (z - 6)   - 108.558739860041 ⎥
⎢                                                                      ⎥
⎢                  ________________________________                    ⎥
⎢                 ╱        2          2          2                     ⎥
⎢ 299792458⋅t + ╲╱  (x - 3)  + (y - 6)  + (z - 9)   - 104.971424683101 ⎥
⎢                                                                      ⎥
⎢                 __________________________________                   ⎥
⎢                ╱        2           2           2                    ⎥
⎣299792458⋅t + ╲╱  (x - 5)  + (y - 10)  + (z - 15)   - 97.8314877736202⎦

In [27]: F.jacobian([x, y, z, t]).subs({y:x, z:x}).det()                                                                          
Out[27]: 0

It is possible to solve this system exactly but a bit of work is needed. We rearrange for the square roots and square everything and then solve (this technique can induce spurious solutions):

In [50]: getsqrt = lambda e: list((a for a in e.atoms(Pow) if a.exp == S.Half))[0]                                                

In [51]: eqs = [getsqrt(e)**2-solve(e, getsqrt(e))[0]**2 for e in F]                                                              

In [52]: solve(eqs, [x, y, z, t])                                                                                                 
Out[52]: 
[(55.9999999982927, 67.0000000004083, 75.9999999981356, 6.05322891533058e-18), (56.6666666664221, 65.6666666641494, 76.666666666
265, 6.05322891533058e-18)]

We can find one of these using nsolve:

In [56]: nsolve(F, [x, y, z, t], [100, 99, 100, 1])                                                                               
Out[56]: 
⎡  56.6666666669817   ⎤
⎢                     ⎥
⎢   65.666666666105   ⎥
⎢                     ⎥
⎢  76.6666666669866   ⎥
⎢                     ⎥
⎣-1.41498771053963e-19⎦

The result for t is different but they are both basically saying that t=0.

In general the reason that this is tricky is that the matrix is badly scaled e.g.:

In [78]: F.jacobian([x, y, z, t]).subs({x:51, y:55, z:60, t:1}).n().eigenvals()                                                   
Out[78]: {-0.0120397251584435: 1, -8.56399289233542e-5: 1, 0.00115898129572641: 1, 299792459.731957: 1}

Newton-Raphson implementation, SymPy linsolve producing empty set?

1 Answers