Python:高斯-牛顿法仅对第一次迭代正确
Python: Gauss-Newton Method only correct for 1st iteration
我正在尝试使用高斯-牛顿法在 Python 中找到 n 个圆的近交点。这里的目标是当给定点 (x1, y1), (x2, y2), ..., (xn, yn) 和半径 R1, R2, . .., Rn,高斯-牛顿法用于找到到n圆的距离平方和最小的点。
这是一种迭代方法,我们从初始向量 v0 = [0, 0] 开始。第一次迭代完全正确,但后续迭代不正确。我找不到代码中的错误:
solutions = []
rmse = []
vk = [[0], [0]] # initial is set here
x = [0, 1, 0]
y = [1, 1, -1]
radii = [1, 1, 1]
iterations = 3
base_str_x = "(x _xi) / sqrt((x xi)^2 + (y yi)^2)"
base_str_y = "(y _yi) / sqrt((x xi)^2 + (y yi)^2)"
base_str_rk = "sqrt((x xi)^2 + (y yi)^2) Ri K"
it = 0
while it < iterations:
i = 0
A = []
while i < len(x):
A.append(["", ""])
A0_str = base_str_x.replace(" xi", "%+f" % (x[i]))
A0_str = A0_str.replace(" yi", "%+f" % (y[i]))
A0_str = A0_str.replace("_xi", "%+f" % (x[i] * -1))
A[i][0] = float(f(vk[0][0], A0_str, vk[1][0]))
A1_str = base_str_y.replace(" xi", "%+f" % (x[i]))
A1_str = A1_str.replace(" yi", "%+f" % (y[i]))
A1_str = A1_str.replace("_yi", "%+f" % (y[i] * -1))
A[i][1] = float(f(vk[0][0], A1_str, vk[1][0]))
i += 1
i = 0
rk = []
while i < len(x):
rk.append([""])
r0_str = base_str_rk.replace(" xi", "%+f" % (x[i]))
r0_str = r0_str.replace(" yi", "%+f" % (y[i]))
r0_str = r0_str.replace("Ri", "%+f" % (radii[i] * -1))
rk[i][0] = float(f(vk[0][0], r0_str, vk[1][0]))
i += 1
lhs = np.matmul(map(list, zip(*A)), A)
rhs = -np.matmul(map(list, zip(*A)), rk)
vk = np.matmul(np.linalg.inv(lhs), rhs)
solutions.append(vk)
此代码计算 A = Dr(x, y),其定义为:
然后求解方程:
对于 vk
。
函数f
是一个"universal function",用于计算A
的每个元素。非常感谢任何有助于找出为什么在第一次迭代之后的后续迭代不正确的帮助。
发生这种情况的原因有两个:
- 括号内的符号是错误的。为了解决这个问题,在每个术语中将
x[1]
乘以 -1
。例如A0_str = A0_str.replace(" yi", "%+f" % (y[i] * -1))
- 错误传播。结果矩阵给出的结果比在后续迭代中转换为
float
的小数位数少。
第一个问题很容易解决,第二个问题就不行了。 matmul
不支持使用 Decimal
并给出 TypeError
。 float64
或 float128
数据类型不出现提高精度。
我正在尝试使用高斯-牛顿法在 Python 中找到 n 个圆的近交点。这里的目标是当给定点 (x1, y1), (x2, y2), ..., (xn, yn) 和半径 R1, R2, . .., Rn,高斯-牛顿法用于找到到n圆的距离平方和最小的点。
这是一种迭代方法,我们从初始向量 v0 = [0, 0] 开始。第一次迭代完全正确,但后续迭代不正确。我找不到代码中的错误:
solutions = []
rmse = []
vk = [[0], [0]] # initial is set here
x = [0, 1, 0]
y = [1, 1, -1]
radii = [1, 1, 1]
iterations = 3
base_str_x = "(x _xi) / sqrt((x xi)^2 + (y yi)^2)"
base_str_y = "(y _yi) / sqrt((x xi)^2 + (y yi)^2)"
base_str_rk = "sqrt((x xi)^2 + (y yi)^2) Ri K"
it = 0
while it < iterations:
i = 0
A = []
while i < len(x):
A.append(["", ""])
A0_str = base_str_x.replace(" xi", "%+f" % (x[i]))
A0_str = A0_str.replace(" yi", "%+f" % (y[i]))
A0_str = A0_str.replace("_xi", "%+f" % (x[i] * -1))
A[i][0] = float(f(vk[0][0], A0_str, vk[1][0]))
A1_str = base_str_y.replace(" xi", "%+f" % (x[i]))
A1_str = A1_str.replace(" yi", "%+f" % (y[i]))
A1_str = A1_str.replace("_yi", "%+f" % (y[i] * -1))
A[i][1] = float(f(vk[0][0], A1_str, vk[1][0]))
i += 1
i = 0
rk = []
while i < len(x):
rk.append([""])
r0_str = base_str_rk.replace(" xi", "%+f" % (x[i]))
r0_str = r0_str.replace(" yi", "%+f" % (y[i]))
r0_str = r0_str.replace("Ri", "%+f" % (radii[i] * -1))
rk[i][0] = float(f(vk[0][0], r0_str, vk[1][0]))
i += 1
lhs = np.matmul(map(list, zip(*A)), A)
rhs = -np.matmul(map(list, zip(*A)), rk)
vk = np.matmul(np.linalg.inv(lhs), rhs)
solutions.append(vk)
此代码计算 A = Dr(x, y),其定义为:
然后求解方程:
对于 vk
。
函数f
是一个"universal function",用于计算A
的每个元素。非常感谢任何有助于找出为什么在第一次迭代之后的后续迭代不正确的帮助。
发生这种情况的原因有两个:
- 括号内的符号是错误的。为了解决这个问题,在每个术语中将
x[1]
乘以-1
。例如A0_str = A0_str.replace(" yi", "%+f" % (y[i] * -1))
- 错误传播。结果矩阵给出的结果比在后续迭代中转换为
float
的小数位数少。
第一个问题很容易解决,第二个问题就不行了。 matmul
不支持使用 Decimal
并给出 TypeError
。 float64
或 float128
数据类型不出现提高精度。