使用 Numpy linalg.lstsq 求解线性系统时获得（明显）不准确的值
Obtaining (significantly) inaccurate values when using Numpy linalg.lstsq to solve a linear system
我正在尝试使用 Numpy 的 linalg.lstsq() 来遵循 this post 中的逻辑。我会在这里重复。
假设我有四个 3 维麦克风 space（或光接收器等）。我想使用每个麦克风记录枪声之间的时间差，或者到达每个麦克风的时间，来确定射手的坐标。所以，我们知道了四个麦克风的坐标x_i, y_i，以及声音的速度v。我们不知道枪手的坐标，也不知道开枪的时间。
对于每个麦克风，我们写：
(X-x_i)**2 + (Y-y_i)**2 + (Z-z_i)**2 = (v(t_i - T))**2
其中X,Y,Z为射手坐标，T为射门时间
如果我们考虑所有 i/j 种可能性，从等式 j 中减去等式 i，问题可以得到简化。我们最终得到 6 方程式（通常，给定 n 个麦克风，n(n-1)/2 方程式）形式为：
2*(x_j - x_i)*X + 2*(y_j - y_i)*Y + 2*(z_j - z_i)*Z + 2 * v**2 * (t_i - t_j) = 2 * v**2 ( t_i**2 - t_j**2) + (x_j**2 + y_j**2 + z_j**2) - (x_i**2 + y_i**2 + z_i**2)
这些方程的形式为 Xv_1 + Y_v2 + Z_v3 + T_v4 = b，其中 v_i 是向量。然后，我们（应该）能够使用线性最小二乘法拟合来“求解”（或者至少粗略估计）X,Y,Z，并且作为副产品，T.
我已经尝试 很多 不同的方法来解决这个问题，所有迹象都指向使用链接 post 中描述的方法（即解决这个问题线性最小二乘法，然后使用此“解决方案”作为更“精确”算法的初始估计）。因此，我编写了以下（不可否认，有些粗糙）Python 代码来开始。它解决了问题的线性最小二乘部分：
from dataclasses import dataclass
from numpy import *
from random import randrange
@dataclass
class Vertexer:
@staticmethod
def find():
# Pick microphones to be at random locations
x_1 = randrange(10); y_1 = randrange(10); z_1 = randrange(10)
x_2 = randrange(10); y_2 = randrange(10); z_2 = randrange(10)
x_3 = randrange(10); y_3 = randrange(10); z_3 = randrange(10)
x_4 = randrange(10); y_4 = randrange(10); z_4 = randrange(10)
# Pick shooter to be at random location
x = randrange(100); y = randrange(100); z = randrange(100)
# Set velocity (ok, it's a ray gun...)
c = 299792 # km/s
# Generate simulated source
t_1 = math.sqrt( (x - x_1)**2 + (y - y_1)**2 + (z - z_1)**2 ) / c
t_2 = math.sqrt( (x - x_2)**2 + (y - y_2)**2 + (z - z_2)**2 ) / c
t_3 = math.sqrt( (x - x_3)**2 + (y - y_3)**2 + (z - z_3)**2 ) / c
t_4 = math.sqrt( (x - x_4)**2 + (y - y_4)**2 + (z - z_4)**2 ) / c
A = array([[2 * (x_1 - x_2), 2 * (y_1 - y_2), 2 * (z_1 - z_2), 2 * c * c * (t_2 - t_1)],
[2 * (x_1 - x_3), 2 * (y_1 - y_3), 2 * (z_1 - z_3), 2 * c * c * (t_3 - t_1)],
[2 * (x_1 - x_4), 2 * (y_1 - y_4), 2 * (z_1 - z_4), 2 * c * c * (t_4 - t_1)],
[2 * (x_2 - x_3), 2 * (y_2 - y_3), 2 * (z_2 - z_3), 2 * c * c * (t_3 - t_2)],
[2 * (x_2 - x_4), 2 * (y_2 - y_4), 2 * (z_2 - z_4), 2 * c * c * (t_4 - t_2)],
[2 * (x_3 - x_4), 2 * (y_3 - y_4), 2 * (z_3 - z_4), 2 * c * c * (t_4 - t_3)]
])
b = array([2 * c * c * (t_2 * t_2 - t_1 * t_1) + (x_1 * x_1 + y_1 * y_1 + z_1 * z_1) - (x_2 * x_2 + y_2 * y_2 + z_2 * z_2),
2 * c * c * (t_3 * t_3 - t_1 * t_1) + (x_1 * x_1 + y_1 * y_1 + z_1 * z_1) - (x_3 * x_3 + y_3 * y_3 + z_3 * z_3),
2 * c * c * (t_4 * t_4 - t_1 * t_1) + (x_1 * x_1 + y_1 * y_1 + z_1 * z_1) - (x_4 * x_4 + y_4 * y_4 + z_4 * z_4),
2 * c * c * (t_3 * t_3 - t_2 * t_2) + (x_2 * x_2 + y_2 * y_2 + z_2 * z_2) - (x_3 * x_3 + y_3 * y_3 + z_3 * z_3),
2 * c * c * (t_4 * t_4 - t_2 * t_2) + (x_2 * x_2 + y_2 * y_2 + z_2 * z_2) - (x_4 * x_4 + y_4 * y_4 + z_4 * z_4),
2 * c * c * (t_4 * t_4 - t_3 * t_3) + (x_3 * x_3 + y_3 * y_3 + z_3 * z_3) - (x_4 * x_4 + y_4 * y_4 + z_4 * z_4)])
solved, residuals, rank, s = linalg.lstsq(A, b)
print(solved)
myVertexer = Vertexer()
myVertexer.find()
不幸的是，它似乎并没有产生非常合理的估计。如果愿意，您可以 运行 自己编写代码，或者这里有一组 运行 将找到的值与算法生成的“模拟”源的坐标进行比较。
Predicted: [ 2.26739519e+00 4.80191502e+00 -5.07181020e+00 ]
Actual: 78 35 57
Predicted: [ 7.72173135e-01 -2.61250803e+00 4.10306750e+00 ]
Actual: 30 75 48
Predicted: [-1.65368110e+00 8.82919123e-01 1.43648336e+00 ]
Actual: 67 56 44
Predicted: [ 4.08698715e+00 -3.19377148e-01 8.81706364e-01 ]
Actual: 9 82 13
诚然，我不确定初始最小二乘解的精度是多少。但是，这些值对我来说似乎不太好。我哪里做错了？
如果我可以提供任何进一步的说明，请告诉我。
Note that I'm aware there are other methods to solving this problem (notably, a pretty good solution can be obtained by simply plugging into one of a number of packages with nonlinear solvers out there). I'm interested in determining why this solution (which I know the linked poster has used successfully, as have several others) is not working, both out of curiosity, and also because I'd like to take this further, where "custom" code will be useful.
有两个不同的问题：
- 你在原文中引用的减法方程有错误post：在RHS上应该是
c * c * (t_i * t_i - t_j * t_j)而不是2 * c * c * (t_i * t_i - t_j * t_j)。
- 第二个问题是在矩阵
A 中有一列 tau （第四列），而在 t_1 到 t_4 的计算中你隐含地假定 tau=0，因此您 不能 将其作为矩阵 A 中的变量；因此需要删除第四列。
修改后的代码如下：
from random import randrange, seed
from numpy import array, linalg, set_printoptions
from dataclasses import dataclass
set_printoptions(suppress=True, linewidth=1000, precision=6)
seed(5)
@dataclass
class Vertexer:
@staticmethod
def find():
# Pick microphones to be at random locations
x_1 = randrange(10); y_1 = randrange(10); z_1 = randrange(10)
x_2 = randrange(10); y_2 = randrange(10); z_2 = randrange(10)
x_3 = randrange(10); y_3 = randrange(10); z_3 = randrange(10)
x_4 = randrange(10); y_4 = randrange(10); z_4 = randrange(10)
# Pick shooter to be at random location
x = randrange(100); y = randrange(100); z = randrange(100)
# Set velocity (ok, it's a ray gun...)
c = 299792 # km/ns
# Generate simulated source
t_1 = math.sqrt( (x - x_1)**2 + (y - y_1)**2 + (z - z_1)**2 ) / c
t_2 = math.sqrt( (x - x_2)**2 + (y - y_2)**2 + (z - z_2)**2 ) / c
t_3 = math.sqrt( (x - x_3)**2 + (y - y_3)**2 + (z - z_3)**2 ) / c
t_4 = math.sqrt( (x - x_4)**2 + (y - y_4)**2 + (z - z_4)**2 ) / c
A = array([[2 * (x_1 - x_2), 2 * (y_1 - y_2), 2 * (z_1 - z_2)],
[2 * (x_1 - x_3), 2 * (y_1 - y_3), 2 * (z_1 - z_3)],
[2 * (x_1 - x_4), 2 * (y_1 - y_4), 2 * (z_1 - z_4)],
[2 * (x_2 - x_3), 2 * (y_2 - y_3), 2 * (z_2 - z_3)],
[2 * (x_2 - x_4), 2 * (y_2 - y_4), 2 * (z_2 - z_4)],
[2 * (x_3 - x_4), 2 * (y_3 - y_4), 2 * (z_3 - z_4)]
])
b = array([1 * c * c * (t_2 * t_2 - t_1 * t_1) + (x_1 * x_1 + y_1 * y_1 + z_1 * z_1) - (x_2 * x_2 + y_2 * y_2 + z_2 * z_2),
1 * c * c * (t_3 * t_3 - t_1 * t_1) + (x_1 * x_1 + y_1 * y_1 + z_1 * z_1) - (x_3 * x_3 + y_3 * y_3 + z_3 * z_3),
1 * c * c * (t_4 * t_4 - t_1 * t_1) + (x_1 * x_1 + y_1 * y_1 + z_1 * z_1) - (x_4 * x_4 + y_4 * y_4 + z_4 * z_4),
1 * c * c * (t_3 * t_3 - t_2 * t_2) + (x_2 * x_2 + y_2 * y_2 + z_2 * z_2) - (x_3 * x_3 + y_3 * y_3 + z_3 * z_3),
1 * c * c * (t_4 * t_4 - t_2 * t_2) + (x_2 * x_2 + y_2 * y_2 + z_2 * z_2) - (x_4 * x_4 + y_4 * y_4 + z_4 * z_4),
1 * c * c * (t_4 * t_4 - t_3 * t_3) + (x_3 * x_3 + y_3 * y_3 + z_3 * z_3) - (x_4 * x_4 + y_4 * y_4 + z_4 * z_4)])
solved, residuals, rank, s = linalg.lstsq(A, b, rcond=None)
print('Actual:', [x, y, z])
print('Predicted:', solved)
myVertexer = Vertexer()
myVertexer.find()
输出为：
Actual: [31, 48, 69]
Predicted: [31. 48. 69.]
我正在尝试使用 Numpy 的 linalg.lstsq() 来遵循 this post 中的逻辑。我会在这里重复。
假设我有四个 3 维麦克风 space（或光接收器等）。我想使用每个麦克风记录枪声之间的时间差，或者到达每个麦克风的时间，来确定射手的坐标。所以，我们知道了四个麦克风的坐标x_i, y_i，以及声音的速度v。我们不知道枪手的坐标，也不知道开枪的时间。
对于每个麦克风，我们写：
(X-x_i)**2 + (Y-y_i)**2 + (Z-z_i)**2 = (v(t_i - T))**2
其中X,Y,Z为射手坐标，T为射门时间
如果我们考虑所有 i/j 种可能性，从等式 j 中减去等式 i，问题可以得到简化。我们最终得到 6 方程式（通常，给定 n 个麦克风，n(n-1)/2 方程式）形式为：
2*(x_j - x_i)*X + 2*(y_j - y_i)*Y + 2*(z_j - z_i)*Z + 2 * v**2 * (t_i - t_j) = 2 * v**2 ( t_i**2 - t_j**2) + (x_j**2 + y_j**2 + z_j**2) - (x_i**2 + y_i**2 + z_i**2)
这些方程的形式为 Xv_1 + Y_v2 + Z_v3 + T_v4 = b，其中 v_i 是向量。然后，我们（应该）能够使用线性最小二乘法拟合来“求解”（或者至少粗略估计）X,Y,Z，并且作为副产品，T.
我已经尝试 很多 不同的方法来解决这个问题，所有迹象都指向使用链接 post 中描述的方法（即解决这个问题线性最小二乘法，然后使用此“解决方案”作为更“精确”算法的初始估计）。因此，我编写了以下（不可否认，有些粗糙）Python 代码来开始。它解决了问题的线性最小二乘部分：
from dataclasses import dataclass
from numpy import *
from random import randrange
@dataclass
class Vertexer:
@staticmethod
def find():
# Pick microphones to be at random locations
x_1 = randrange(10); y_1 = randrange(10); z_1 = randrange(10)
x_2 = randrange(10); y_2 = randrange(10); z_2 = randrange(10)
x_3 = randrange(10); y_3 = randrange(10); z_3 = randrange(10)
x_4 = randrange(10); y_4 = randrange(10); z_4 = randrange(10)
# Pick shooter to be at random location
x = randrange(100); y = randrange(100); z = randrange(100)
# Set velocity (ok, it's a ray gun...)
c = 299792 # km/s
# Generate simulated source
t_1 = math.sqrt( (x - x_1)**2 + (y - y_1)**2 + (z - z_1)**2 ) / c
t_2 = math.sqrt( (x - x_2)**2 + (y - y_2)**2 + (z - z_2)**2 ) / c
t_3 = math.sqrt( (x - x_3)**2 + (y - y_3)**2 + (z - z_3)**2 ) / c
t_4 = math.sqrt( (x - x_4)**2 + (y - y_4)**2 + (z - z_4)**2 ) / c
A = array([[2 * (x_1 - x_2), 2 * (y_1 - y_2), 2 * (z_1 - z_2), 2 * c * c * (t_2 - t_1)],
[2 * (x_1 - x_3), 2 * (y_1 - y_3), 2 * (z_1 - z_3), 2 * c * c * (t_3 - t_1)],
[2 * (x_1 - x_4), 2 * (y_1 - y_4), 2 * (z_1 - z_4), 2 * c * c * (t_4 - t_1)],
[2 * (x_2 - x_3), 2 * (y_2 - y_3), 2 * (z_2 - z_3), 2 * c * c * (t_3 - t_2)],
[2 * (x_2 - x_4), 2 * (y_2 - y_4), 2 * (z_2 - z_4), 2 * c * c * (t_4 - t_2)],
[2 * (x_3 - x_4), 2 * (y_3 - y_4), 2 * (z_3 - z_4), 2 * c * c * (t_4 - t_3)]
])
b = array([2 * c * c * (t_2 * t_2 - t_1 * t_1) + (x_1 * x_1 + y_1 * y_1 + z_1 * z_1) - (x_2 * x_2 + y_2 * y_2 + z_2 * z_2),
2 * c * c * (t_3 * t_3 - t_1 * t_1) + (x_1 * x_1 + y_1 * y_1 + z_1 * z_1) - (x_3 * x_3 + y_3 * y_3 + z_3 * z_3),
2 * c * c * (t_4 * t_4 - t_1 * t_1) + (x_1 * x_1 + y_1 * y_1 + z_1 * z_1) - (x_4 * x_4 + y_4 * y_4 + z_4 * z_4),
2 * c * c * (t_3 * t_3 - t_2 * t_2) + (x_2 * x_2 + y_2 * y_2 + z_2 * z_2) - (x_3 * x_3 + y_3 * y_3 + z_3 * z_3),
2 * c * c * (t_4 * t_4 - t_2 * t_2) + (x_2 * x_2 + y_2 * y_2 + z_2 * z_2) - (x_4 * x_4 + y_4 * y_4 + z_4 * z_4),
2 * c * c * (t_4 * t_4 - t_3 * t_3) + (x_3 * x_3 + y_3 * y_3 + z_3 * z_3) - (x_4 * x_4 + y_4 * y_4 + z_4 * z_4)])
solved, residuals, rank, s = linalg.lstsq(A, b)
print(solved)
myVertexer = Vertexer()
myVertexer.find()
不幸的是，它似乎并没有产生非常合理的估计。如果愿意，您可以 运行 自己编写代码，或者这里有一组 运行 将找到的值与算法生成的“模拟”源的坐标进行比较。
Predicted: [ 2.26739519e+00 4.80191502e+00 -5.07181020e+00 ]
Actual: 78 35 57
Predicted: [ 7.72173135e-01 -2.61250803e+00 4.10306750e+00 ]
Actual: 30 75 48
Predicted: [-1.65368110e+00 8.82919123e-01 1.43648336e+00 ]
Actual: 67 56 44
Predicted: [ 4.08698715e+00 -3.19377148e-01 8.81706364e-01 ]
Actual: 9 82 13
诚然，我不确定初始最小二乘解的精度是多少。但是，这些值对我来说似乎不太好。我哪里做错了？
如果我可以提供任何进一步的说明，请告诉我。
Note that I'm aware there are other methods to solving this problem (notably, a pretty good solution can be obtained by simply plugging into one of a number of packages with nonlinear solvers out there). I'm interested in determining why this solution (which I know the linked poster has used successfully, as have several others) is not working, both out of curiosity, and also because I'd like to take this further, where "custom" code will be useful.
有两个不同的问题：
- 你在原文中引用的减法方程有错误post：在RHS上应该是
c * c * (t_i * t_i - t_j * t_j)而不是2 * c * c * (t_i * t_i - t_j * t_j)。 - 第二个问题是在矩阵
A中有一列tau（第四列），而在t_1到t_4的计算中你隐含地假定tau=0，因此您 不能 将其作为矩阵A中的变量；因此需要删除第四列。
修改后的代码如下：
from random import randrange, seed
from numpy import array, linalg, set_printoptions
from dataclasses import dataclass
set_printoptions(suppress=True, linewidth=1000, precision=6)
seed(5)
@dataclass
class Vertexer:
@staticmethod
def find():
# Pick microphones to be at random locations
x_1 = randrange(10); y_1 = randrange(10); z_1 = randrange(10)
x_2 = randrange(10); y_2 = randrange(10); z_2 = randrange(10)
x_3 = randrange(10); y_3 = randrange(10); z_3 = randrange(10)
x_4 = randrange(10); y_4 = randrange(10); z_4 = randrange(10)
# Pick shooter to be at random location
x = randrange(100); y = randrange(100); z = randrange(100)
# Set velocity (ok, it's a ray gun...)
c = 299792 # km/ns
# Generate simulated source
t_1 = math.sqrt( (x - x_1)**2 + (y - y_1)**2 + (z - z_1)**2 ) / c
t_2 = math.sqrt( (x - x_2)**2 + (y - y_2)**2 + (z - z_2)**2 ) / c
t_3 = math.sqrt( (x - x_3)**2 + (y - y_3)**2 + (z - z_3)**2 ) / c
t_4 = math.sqrt( (x - x_4)**2 + (y - y_4)**2 + (z - z_4)**2 ) / c
A = array([[2 * (x_1 - x_2), 2 * (y_1 - y_2), 2 * (z_1 - z_2)],
[2 * (x_1 - x_3), 2 * (y_1 - y_3), 2 * (z_1 - z_3)],
[2 * (x_1 - x_4), 2 * (y_1 - y_4), 2 * (z_1 - z_4)],
[2 * (x_2 - x_3), 2 * (y_2 - y_3), 2 * (z_2 - z_3)],
[2 * (x_2 - x_4), 2 * (y_2 - y_4), 2 * (z_2 - z_4)],
[2 * (x_3 - x_4), 2 * (y_3 - y_4), 2 * (z_3 - z_4)]
])
b = array([1 * c * c * (t_2 * t_2 - t_1 * t_1) + (x_1 * x_1 + y_1 * y_1 + z_1 * z_1) - (x_2 * x_2 + y_2 * y_2 + z_2 * z_2),
1 * c * c * (t_3 * t_3 - t_1 * t_1) + (x_1 * x_1 + y_1 * y_1 + z_1 * z_1) - (x_3 * x_3 + y_3 * y_3 + z_3 * z_3),
1 * c * c * (t_4 * t_4 - t_1 * t_1) + (x_1 * x_1 + y_1 * y_1 + z_1 * z_1) - (x_4 * x_4 + y_4 * y_4 + z_4 * z_4),
1 * c * c * (t_3 * t_3 - t_2 * t_2) + (x_2 * x_2 + y_2 * y_2 + z_2 * z_2) - (x_3 * x_3 + y_3 * y_3 + z_3 * z_3),
1 * c * c * (t_4 * t_4 - t_2 * t_2) + (x_2 * x_2 + y_2 * y_2 + z_2 * z_2) - (x_4 * x_4 + y_4 * y_4 + z_4 * z_4),
1 * c * c * (t_4 * t_4 - t_3 * t_3) + (x_3 * x_3 + y_3 * y_3 + z_3 * z_3) - (x_4 * x_4 + y_4 * y_4 + z_4 * z_4)])
solved, residuals, rank, s = linalg.lstsq(A, b, rcond=None)
print('Actual:', [x, y, z])
print('Predicted:', solved)
myVertexer = Vertexer()
myVertexer.find()
输出为：
Actual: [31, 48, 69]
Predicted: [31. 48. 69.]