你如何处理 Matlab 中的精度问题？

Question

我正在尝试在 Matlab 中编写一个生日问题计算器，但在 (1 - 非常小的浮点数 = 1) 处存在精度问题。

我目前的问题是，我想看看在有 23,000,000 个活动会话令牌的网站上猜测 UUID 需要多少次尝试，这些令牌有 128 位可能的唯一值，因此猜测有效会话令牌的几率令牌超过 50%。

我首先通过以下方式模拟该过程：

1. 我将 success_rate 设置为 (23,000,000 / (2^128))
2. 我将 failure_rate 设置为 (1 - success_rate)

但后来我注意到这个值为 1。

更糟糕的是，将(1 - 23000000/(2^128))^n > 0.5 输入Wolfram Alpha 并没有提供有用的答案。

我的第一个想法是完全抛弃 Matlab 并在 Java 中创建自己的库，它根本不使用浮点值，而是将比率存储为 BigDecimal 对象对，这将通过仅进行计算来消除所有精度问题在最后一点，并将此计算存储为一对最小值 - 最大值，以将结果显示为解决方案所在的范围（其中不存在精确解决方案，因为浮点除法会导致错误和无法使用表示的值指定精度的浮点数，但可以通过仅指定精确的实际比率来表示准确的答案，因为从不对其应用除法，而是显示比率）。

有没有办法在不必发明这样的系统的情况下解决这类问题，或者这些问题本质上是不可能使用浮点系统解决的？

Answer 1

...这些问题是否天生就无法使用浮点系统解决？

简短说明：

嗯，在 MATLAB 中默认是可以的，如果您在 MATLAB 中使用符号工具箱，则不是。

您绝对可以在 MATLAB 中用双精度浮点数表示非常小的数字。但是，您遇到的问题与对彼此相差太多数量级的双精度浮点数进行运算有关 - 在执行计算时，您会受到 MATLAB 计算精度的限制。

谢天谢地，有一个以符号工具箱和variable-precision arithmetic 的形式缓解此问题的工具箱。如果您想在执行1 - (small_value) 时得到 1 以外的值，请查看它。

更长的解释：

http://www.mathworks.com/help/matlab/matlab_prog/floating-point-numbers.html#f2-98720

MATLAB 中的双精度浮点数具有令人印象深刻的最大精度 -1.79769e+308 to -2.22507e-308 and 2.22507e-308 to 1.79769e+308。但是，MATLAB 最多只能计算 53 位精度：精度为 9.007199255×10¹⁵。

这是我对如何产生您遇到的结果的解释 (1 - small_value = 1)：

数字1.234e12 以大约1e16 的精度表示，这意味着MATLAB 可以对这个数字进行操作，误差大约为1e-4。同样，2.345e-7 的计算误差大致为1e-23。因此，将这两个数字相加将产生1e-4 的错误，因此较小的数字已在 MATLAB 执行的计算错误中丢失。

如果您不介意等待更长的计算时间来执行比 53 位大得多的操作，那么我强烈建议您使用 MATLAB 中的符号工具箱（即vpa 函数）。

如果我的回答不适合您，也许您可​​以在 MATLAB 论坛中查看此 answer to a related question。我从这个答案中提取了部分样本编号。

编码愉快，希望对您有所帮助！

Answer 2

简单解释：

使用：

   eps(double(1))

在 Matlab 中，您会发现 1（最大精度 = 双倍）和它在执行数学运算时可以区分的下一个浮点数之间的最小差距。在这种情况下，差距等于 2.2204e-016

自：

success_rate = (23,000,000 / (2^128))

将返回6.7591e-032，并且在执行 1 - 6.7591e-032 时它比上面引入的间隙要小得多，Matlab 理解这是从 1 中减去 0，因此你总是得到 1 作为答案。希望对您有所帮助。

Answer 3

正如所有其他答案所指出的那样，问题是r = 3000000/(2^128) < eps(1)/2，所以1 + r == 1

最简单的方法是重新排列表达式，并在此过程中利用其他一些功能。重写：

(1 - 23000000/(2^128))^n = exp(n*log(1- 23000000/(2^128))

现在，这仍然会出现同样的问题，但是有一个log1p 函数可以准确计算log(1+x)。所以改为使用：

exp(n*log1p(-23000000/(2^128)))

Answer 4

其他答案已经解释了为什么您无法执行所需的计算，因为您使用的数字大小存在差异。但是，正如我在评论中提到的，您可以尝试使用较小的数字来显示趋势。我们将“预计”值称为size_of_key_space / (2 * number_of_keys)。对于获得 50% 的成功概率，这是一种幼稚的期望值。为了证明这一点，我对许多不同的键集和键空间进行了模拟。都很大，具有不同的稀疏性：

function sparse_probability()

num_keys = logspace(2, 5, 15);  % number of keys varies from 1e2 to 1e5
key_spaces = logspace(6, 12, 15);  % size of key space varies from 1e6 to 1e12
% so p_sucess varies from 1e-4 to 1e-7

num_experiments = length(num_keys);

results = zeros(1,num_experiments);
proportions = zeros(1,num_experiments);

for i = 1:num_experiments
    num_objs = num_keys(i);
    size_of_key_space = key_spaces(i);
    p_success = num_objs/size_of_key_space;
    p_fail = 1 - p_success;

    total_fail = 1;
    num_trials = 0;
    while (total_fail > 0.5)
        total_fail = total_fail * p_fail;
        num_trials = num_trials + 1;
    end


    results(i) = num_trials;
    proportions(i) = num_trials/(size_of_key_space/(2*num_objs));
    fprintf('p_success = %f, num_trials = %d, ratio = %f, num_keys = %e; size key_space = %e\n', 1 - total_fail, num_trials, proportions(i), num_objs, size_of_key_space);
end

由于键集和键空间的大小差异很大，我计算了上面“预计”值的比率，以及达到 50% 概率所需的实际试验次数。上面函数的输出是：

p_success = 0.500044, num_trials = 6932, ratio = 1.386400, num_keys = 1.000000e+02; size key_space = 1.000000e+06
p_success = 0.500010, num_trials = 11353, ratio = 1.386293, num_keys = 1.637894e+02; size key_space = 2.682696e+06
p_success = 0.500006, num_trials = 18595, ratio = 1.386292, num_keys = 2.682696e+02; size key_space = 7.196857e+06
p_success = 0.500008, num_trials = 30457, ratio = 1.386309, num_keys = 4.393971e+02; size key_space = 1.930698e+07
p_success = 0.500004, num_trials = 49885, ratio = 1.386300, num_keys = 7.196857e+02; size key_space = 5.179475e+07
p_success = 0.500001, num_trials = 81706, ratio = 1.386294, num_keys = 1.178769e+03; size key_space = 1.389495e+08
p_success = 0.500001, num_trials = 133826, ratio = 1.386297, num_keys = 1.930698e+03; size key_space = 3.727594e+08
p_success = 0.500002, num_trials = 219193, ratio = 1.386298, num_keys = 3.162278e+03; size key_space = 1.000000e+09
p_success = 0.500001, num_trials = 359014, ratio = 1.386295, num_keys = 5.179475e+03; size key_space = 2.682696e+09
p_success = 0.500001, num_trials = 588027, ratio = 1.386296, num_keys = 8.483429e+03; size key_space = 7.196857e+09
p_success = 0.500000, num_trials = 963125, ratio = 1.386295, num_keys = 1.389495e+04; size key_space = 1.930698e+10
p_success = 0.500000, num_trials = 1577496, ratio = 1.386294, num_keys = 2.275846e+04; size key_space = 5.179475e+10
p_success = 0.500000, num_trials = 2583771, ratio = 1.386294, num_keys = 3.727594e+04; size key_space = 1.389495e+11
p_success = 0.500000, num_trials = 4231943, ratio = 1.386295, num_keys = 6.105402e+04; size key_space = 3.727594e+11
p_success = 0.500000, num_trials = 6931472, ratio = 1.386294, num_keys = 1.000000e+05; size key_space = 1.000000e+12

如果您要绘制比率列与键空间大小的关系，您会得到一条直线。例如，只要密钥集和密钥空间相隔几个数量级，该比率基本上是恒定的。请注意，稀疏度会有所不同，但这不会影响比率。这是这类稀疏概率问题的典型特征。因此，从这个简单的实验中，您可以非常自信地说，2.3e7 键所需的猜测次数，在2^128 = 3.4e38 的键空间中，是高于1.386294 的比率限制与预测值的乘积一共

1.386294 * (2^128 / (2 * 2.3e7)) = 1.02550305123542e+31 

猜测有效 UUID 的概率为 50%。

以每秒 1 万亿次的猜测，进行这么多猜测需要 3250 亿年。换句话说，你是安全的。 :)

Answer 5

正如其他人所解释的，(1 - 23000000/2^128) 太接近于 1，无法以双精度浮点值的尾数的 53 位表示，所以 (1 - 230000000/2^128)^ n 无法计算。

其他软件包（python+sympy,mathematica,...）可以进行任意精度计算，matlab有一个多精度计算工具箱。这将允许您直接执行计算。

您可以将方程重新排列为二项式展开：

(a + b)^n = a^n + C(1,n)a^(n-1)b + C(2,n)a^(n-2)b^2 + ...

其中 C(k,n) 是从大小为 n 的池中选择 k 个项目的方法数。由于b^k 对于较大的k 来说很小，因此请忽略这些术语，并将其近似为：

(1 - b)^n = 1 - n b + O(b^2)

b = 23000000/2^38。为n 求解1 - n b = 0.5 得到其他人给出的近似n = 2^128 / (2 * 23000000)。

Herbie 有时可以帮助您重写方程以提高数值稳定性。

另一个最受欢迎的技巧是在您试图逼近的值附近执行泰勒展开，给出一个可以在一系列输入中使用的多项式。可以使用多精度库确定多项式次数和有效范围，以便您知道您的值在整个范围内都准确到机器精度。 Wolfram Alpha 提供了一个在线泰勒级数计算器。

更多细节可以在以下书籍中找到：

1. 新泽西州海厄姆。数值算法的准确性和稳定性：第二版。暹; 2002 年。