可靠地确定浮点值向量是不是（数字上）等距

Posted 2023-04-14

【中文标题】可靠地确定浮点值向量是不是（数字上）等距

【英文标题】：Reliably determine if vector of floating point values is (numerically) equidistant可靠地确定浮点值向量是否（数字上）等距

【发布时间】：2020-07-15 13:47:06

【问题描述】：

假设有一个浮点值向量。

X = [x1, x2, x3 ... xN]

我想检查这些值是否等距分布。当然，由于浮点运算，不能保证相邻的差值都是一样的。

执行此检查的稳健算法是什么？

从某种意义上说是强大的

不对向量元素的大小和符号做任何假设 也不是向量的长度 无论等距序列是通过重复递增还是通过乘法i*dx + x0 生成的。

目前我查看差异向量，将它们四舍五入到一定数量的小数位，然后检查它们是否都相同。然而，这不够健壮，如果向量很长，有时会失败。 我选择的小数位数是最大绝对向量值的log10乘以eps乘以10，即log10( max(abs(vec))*eps*10 )

编辑：这是当前代码（在 MATLAB 中）和失败的情况。一切都是doublei.e. IEEE binary64 类型

function [b, v0, dv] = isequidistant(vec)
    diffv = diff(vec);
    % Note round(x,N) rounds to the N-th decimal place, e.g. round(1234, -2) == 1200
    diffv = round( diffv, -ceil(log10(max(abs(vec(~isinf(vec)&vec ~= 0)))*eps*10)) );
    diffv = unique( diffv );
    if length(diffv) == 1
        b = true;
        v0 = vec(1);
        dv = diffv;
    else
        b = false;
        v0 = vec(1);
        dv = nan;
    end
end

这里失败了： isequidistant( [0:330000]*3.027872455e-01 ) == 0

最有趣的是 isequidistant( [0:330000]*3.02787245e-01 ) == 1

我知道人们可能会玩弄常数和因子，但我在这里寻求更好的解决方案。上面的那个会产生一个虚假的失败，这让我很担心。

【问题讨论】：

你能举一个你的方法失败的输入例子吗（例如，给出一个生成这样一个序列的算法）？重复增量和i*dx+x0 和(x0*(N-i) + xN*i) / N（以及其他替代方法）可能会生成无法通过其他任何一种方法复制的序列。您是否假设 IEEE 754 binary64/binary32 浮点数或任何值（即使是任意基数？）

Re“我选择的小数位数是最大绝对向量值乘以 eps 乘以 10 的 log10”：什么？因此，如果 100 是最大绝对向量值，则 log10 为 2，对于 IEEE-754 binary64，epsilon 为 2^-52，而 10，因此 2•2^-52•10 = 4.44e-15 个小数位。这是没有意义的。再试一次？

你需要的最大 N 是多少？你知道 x0 和 dx 的界限是多少？它们总是非负的吗？任何接受等间距的点但计算中的舍入误差的测试必须接受一些不等间距的点。你对这样的误报有多大的容忍度？什么会导致您的点不等间距 - 它们是否来自某个过程，可以为我们提供关于它们必须偏离等间距多少的信息？

Re “......并且不假设这个等距序列是如何生成的”：这是一个不可能的条件。通常，算术误差的复合会产生范围从 0 到无穷大的误差，或者可以是 NaN。您必须对数字的产生方式有一些限制，以便确定潜在错误的界限。

@chtz 我添加了我的 (MATLAB) 代码和一个失败的测试用例。

【参考方案1】：

一种直接的方法是在第一个和最后一个元素之间生成一个线性间隔序列，并将每个元素与之进行比较。

下面的代码允许与一些容差进行比较。容差参数可以设为可选（例如，默认为0 或eps）：

function [b, v0, dv] = isequidistant(vec, tolerance)
    v0 = vec(1); vN = vec(end); N = length(vec);
    b = all(abs(linspace(v0,vN,N)-vec) <= tolerance * max(abs(v0),abs(vN)));
    dv = (vN - v0) / (N-1);
end

【参考方案2】：

由于允许存在一定程度的不精确性并且仍然认为“等距”，请考虑计算所有差值，记录最高和最低值，然后最后获取 hi 和 lo 相对于平均值的差值。

// Pseudo code
X = [x1, x2, x3 ... xN]
if (N <= 2) return true;
avg = sum(x1...xN)/N
lo = x2 - x1
hi = x2 - x1
for (i = 3 to N)
   diff = xi - x(i-1)
   if (diff < lo) lo = diff
   if (diff > hi) hi = diff
 
 Now assess if lo,hi are too far apart
 Maybe by comparing the |(hi-lo)/avg| <= tolerance