Newton-Raphson 实现，SymPy linsolve 产生空集？答案

【问题标题】：Newton-Raphson implementation, SymPy linsolve producing empty set?Newton-Raphson 实现，SymPy linsolve 产生空集？
【发布时间】：2020-11-27 16:32:44
【问题描述】：

我写了一个牛顿-拉夫森方法（牛顿方法的多元推广）的实现来解决这个系统：

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)

未知数是x, y, z 和t。该系统描述了一个定位问题的解决方案，其中三个观察者的位置、某个信号的速度（在本例中为c = 299792458）以及信号到达每个观察者的时间是已知的，我们想要确定信号源的坐标x, y, z（在此过程中，我们还找到了发射时间t）。

我相信能很好地理解算法。在这个解决系统的小例子中：

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

即使最初的猜测不太好，我也获得了出色的收敛性。

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()

我的代码：

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()

我用一些随机坐标“模拟”了一个源，并选择了观察者的坐标为：

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

产生到达时间：

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

我的问题：

上面的代码抛出错误：

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

似乎linsolve 正在生成EmptySet，但我不确定为什么？请注意，这不应该是该系统的可解性问题。据我了解，Newton-Raphson 是用于解决此问题的“标准”技术。

编辑：

看来A是：

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]])

而b 是：

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

所以A 是单数...不知道我做错了什么？

更新：

我取得了一些的进步。这些时候我喂了算法：

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

对应的坐标：

83 81 76

也就是说，这就是我希望我的 NR 代码生成的内容。通过将我的初始值更改为100, 100, 100（我认为这是相当好的起点），我大部分时间都能够避免奇异矩阵问题。然而，现在，我离我期望的那种价值观还差得很远。同样，我想要83 81 76，但我得到了（每一行都是牛顿方法中的一个步骤）：

(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

这……不好。

【问题讨论】：

标签： python algorithm sympy

【解决方案1】：

如果 x、y 和 z 都相等，则矩阵是奇异的：

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

完全可以解决这个系统，但需要做一些工作。我们重新排列平方根并对所有内容进行平方，然后求解（这种技术会导致虚假解）：

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)]

我们可以使用nsolve找到其中之一：

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

t 的结果不同，但它们基本上都在说t=0。

一般来说，这很棘手的原因是矩阵缩放不当，例如：

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}

【讨论】：