关于几何中的数值稳定性，这是我在Stack Exchange上的第一个问题。我编写了下面这个函数来计算从一个向量到另一个向量的有符号角度，灵感来自于这里提出的想法，我对此深感感激。
在这个函数中，三重积用于从第一部分得到的原始锐角转化为从第一个向量到第二个向量打开的有向角度。
当在上面的示例测试代码中，给定的第二个向量的角度非常接近于第一个给定向量的角度，无论是从其两侧，而不是得到结果角度分别非常接近于0和非常接近于2，函数的结果实际上是零，因此这两个语义上不同的向量方向产生相同的输出角度。
在这种意义上，您将如何使该函数在数值上稳定？
假设有向角度从0增长到接近2的数字，您将如何避免函数的输出角度在无限接近完整圆形有向角度时跳至0？只要其输出是无限精确的，我们就没问题，但如果它从接近2变为0，当它无限接近2时――这种不连续的值跳跃对于我的应用程序以及我猜其他许多情况下对其输出的解释都是灾难性的。
显然，我们可以通过实证找到接近2的角度函数输出为0的上界边界――实际上下面附有一个这样的函数――但是当输入是产生该角度的向量时，这一点毫无帮助（我们无法从向量中得知它将产生一个翻转到0的角度，因此当函数的结果角度为0时――我们无法确定真实的角度是0还是接近2的值）。
我想你可以用任何一个词来替换“有向角”的术语，比如“顺时针（或逆时针）角度”。
以下是诊断函数。
而且说了这么多，我也很高兴能够了解可能在这个函数中出现的其他数值稳定性问题。我看到有关更喜欢使用反正切函数而不是反余弦函数及其变体的评论，但不知道它们是否完全适用于numpy，是否会引入其他稳定性缺陷。
同时，非常感谢对任何关于思考数值稳定性的一般方法论或特定于numpy的安全编写数值稳定计算的指导方针的评论；尽管我意识到不同的应用程序将担心不同的范围，在这些范围内的不准确性对它们的特定逻辑产生不利影响。
由于我输入到这个函数中的数据来自机器学习预测，即使在单次调用函数时它们的概率很低，我也无法避免极端向量值的出现。
谢谢！
在这个函数中，三重积用于从第一部分得到的原始锐角转化为从第一个向量到第二个向量打开的有向角度。
import numpy as np
from numpy import ndarray, sign
from math import pi
Vec3D = ndarray
def directed_angle(
from_vec: Vec3D,
to_vec: Vec3D,
normal: Vec3D,
debug=False):
""" returns the [0, 2) directed angle opening up from the first to the second vector,
given a 3rd vector linearly independent of the first two vectors, which is used to provide
the directionality of the acute (regular, i.e. undirected) angle firstly computed """
assert np.any(from_vec) and np.any(to_vec), \
f'from_vec, to_vec, and normal must not be zero vectors, but from_vec: {from_vec}, to_vec: {to_vec}'
magnitudes = np.linalg.norm(from_vec) * np.linalg.norm(to_vec)
dot_product = np.matmul(from_vec, to_vec)
angle = np.arccos(np.clip(dot_product / magnitudes, -1, 1))
# the triplet product reflects the sign of the angle opening up from the first vector to the second,
# based on the input normal value providing the directionality which is otherwise totally abscent
triple_product = np.linalg.det(np.array([from_vec, to_vec, normal]))
if triple_product < 0:
result_angle = angle
else:
result_angle = 2*pi - angle
# flip the output range from (0, 2] to [0, 2) which is the prefernce for our semantics
if result_angle == 2*pi:
result_angle = 0
if debug:
print(f'\nfrom_vec: {from_vec}, '
f'\nto_vec: {to_vec}'
f'\nnormal: {normal}, '
f'\nundirected angle: {angle},'
f'\nundirected angle/pi: {angle/pi}, '
f'\ntriple product: {triple_product}'
f'\ntriple product sign: {sign(triple_product)}'
f'\nresult_angle/pi: {result_angle/pi}')
return result_angle
def test():
v0 = np.array([1, 0, 0])
normal = np.array([0, 0, 1])
vecs = dict(
same_direction = np.array([1, 0, 0]),
opposite_direction = np.array([-1, 0, 0]),
close_to_same_direction_from_side_1 = np.array([1, 0.00000001, 0]),
close_to_same_direction_from_side_2 = np.array([1, -0.00000001, 0]))
for desc, v in vecs.items():
print(f'\n{desc}:')
directed_angle(v0, v, normal, debug=True)
当在上面的示例测试代码中，给定的第二个向量的角度非常接近于第一个给定向量的角度，无论是从其两侧，而不是得到结果角度分别非常接近于0和非常接近于2，函数的结果实际上是零，因此这两个语义上不同的向量方向产生相同的输出角度。
在这种意义上，您将如何使该函数在数值上稳定？
假设有向角度从0增长到接近2的数字，您将如何避免函数的输出角度在无限接近完整圆形有向角度时跳至0？只要其输出是无限精确的，我们就没问题，但如果它从接近2变为0，当它无限接近2时――这种不连续的值跳跃对于我的应用程序以及我猜其他许多情况下对其输出的解释都是灾难性的。
显然，我们可以通过实证找到接近2的角度函数输出为0的上界边界――实际上下面附有一个这样的函数――但是当输入是产生该角度的向量时，这一点毫无帮助（我们无法从向量中得知它将产生一个翻转到0的角度，因此当函数的结果角度为0时――我们无法确定真实的角度是0还是接近2的值）。
我想你可以用任何一个词来替换“有向角”的术语，比如“顺时针（或逆时针）角度”。
以下是诊断函数。
def test_monotonicity(debug=False, interpolation_interval=1e-4):
""" test that the function is monotonic (up to the used interpolation interval),
and correct, across the interval of its domain ([0, 2)) """
v0 = np.array([1, 0, 0])
normal = np.array([0, 0, 1])
# range from 0 to 2 at 0.1 intervals
angles = np.arange(0, 2*pi, step=interpolation_interval)
prev_angle = None
for angle in angles:
# make a vector representing the current angle away from v0, such that it goes clockwise
# away from v0 as `angle` increases (clockwise assuming you look at the two vectors from
# the perspective of the normal vector, otherwise directionality is void)
vec = np.array([np.cos(angle), -np.sin(angle), 0])
result_angle = directed_angle(v0, vec, normal)
if debug: print(f'angle/pi: {angle/pi}, computed angle/pi: {result_angle/pi}')
if prev_angle is None:
assert angle == 0
else:
assert angle > prev_angle
drift = result_angle - angle
if angle == result_angle:
pass
else:
print(f'angle/pi: {angle/pi}, result_angle/pi: {result_angle/pi}, drift/pi: {drift/pi}')
assert isclose(result_angle, angle, rel_tol=1e-6)
prev_angle = angle
def test_demonstrating_the_concern():
""" demonstration of sufficiently small angles from both sides of the reference vector (v0)
collapsing to zero, which is a problem only when the angle should be close to 2 pi """
v0 = np.array([1, 0, 0])
normal = np.array([0, 0, 1])
vecs = dict(
same_direction = np.array([1, 0, 0]),
opposite_direction = np.array([-1, 0, 0]),
close_to_same_direction_from_side_1 = np.array([1, +0.00000001, 0]),
close_to_same_direction_from_side_2 = np.array([1, -0.00000001, 0]))
expected_results = [0, 1*pi, None, None]
for (desc, v), expected_result in zip(vecs.items(), expected_results):
print(f'\n{desc}:')
v = v / np.linalg.norm(v)
result = directed_angle(v0, v, normal, debug=True)
if expected_result is not None:
assert(result == expected_result)
def test_angle_approaching_2pi(epsilon=1e-7, interpolation_interval=1e-10, verbose=True):
""" stress test the angle values approaching 2pi where the result flips to zero """
v0 = np.array([1, 0, 0])
normal = np.array([0, 0, 1])
# range from 2-epsilon to close to 2
angles = np.arange(2*pi-epsilon, 2*pi, step=interpolation_interval)
prev_angle = None
for angle in angles:
# vector representing the current angle away from v0, clockwise
vec = np.array([np.cos(angle), -np.sin(angle), 0])
result_angle = directed_angle(v0, vec, normal)
if prev_angle is None:
pass
else:
assert angle > prev_angle
drift = result_angle - angle
if angle == result_angle:
pass
else:
if verbose: print(f'angle/pi: {angle / pi}, result_angle/pi: {result_angle / pi}, drift/pi: {drift / pi}')
if result_angle == 0:
print(f'function angle output hit 0 at angle: {angle};'
f'\nangle/pi: {angle / pi}'
f'\nangle distance from 2: {2*pi - angle}')
break
else:
assert isclose(result_angle, angle, rel_tol=1e-6)
prev_angle = angle
而且说了这么多，我也很高兴能够了解可能在这个函数中出现的其他数值稳定性问题。我看到有关更喜欢使用反正切函数而不是反余弦函数及其变体的评论，但不知道它们是否完全适用于numpy，是否会引入其他稳定性缺陷。
同时，非常感谢对任何关于思考数值稳定性的一般方法论或特定于numpy的安全编写数值稳定计算的指导方针的评论；尽管我意识到不同的应用程序将担心不同的范围，在这些范围内的不准确性对它们的特定逻辑产生不利影响。
由于我输入到这个函数中的数据来自机器学习预测，即使在单次调用函数时它们的概率很低，我也无法避免极端向量值的出现。
谢谢！
np.linalg.norm进行归一化步骤：np.linalg.norm(np.array([1, 1e-8, 0],)) == np.linalg.norm(np.array([1, 0, 0]))结果为True。如果是1e-7或者"更小"，则结果为False。 - cardsnp.linalg.norm：np.linalg.norm(np.array([1, 1e-8, 0],)) == np.linalg.norm(np.array([1, 0, 0]))结果为True。如果是1e-7或 "更小" 则结果为False。 - cardsnp.linalg.norm：np.linalg.norm(np.array([1, 1e-8, 0],)) == np.linalg.norm(np.array([1, 0, 0]))结果为True。如果是1e-7或者更小，则结果为False。 - undefined