1

当我试图在 matlab 中模拟我的正弦近似时，我发现了一个奇怪的问题。问题是，当我将函数应用于数组时，它会返回一个结果，而将函数应用于单个值会给出稍微不同的结果。

在此示例中，我能够获得相同的行为：

z = single(0:0.001:1);
F = @(x) (x.^2 - single(1.2342320e-001)).*x.^2;  %some test function

z(999)        % Returns 9.9800003e-001
F(z(999))     % Returns 8.6909407e-001
temp = F(z);  
temp(999)     % Voila! It returns 8.6909401e-001

我也发现了一些东西。一个是第一个结果是正确的（不是后者）。其次，术语的重新排列有时可以解决问题。所以我不知道如何摆脱它。

4

5 回答 5

3

单精度只有最多 7 位十进制数字的有意义的精度，所以说8.6909407e-001“正确”和8.6909401e-001“错误”在这种情况下没有太大意义。

正如您已经发现的那样，浮点算术对运算顺序也很敏感。当对矩阵而不是标量进行操作时，Matlab 可能会巧妙地改变计算顺序。

于 2012-06-09T19:22:16.073 回答
1

对于向量算术，MatLab 的线性代数库可能使用 SIMD 指令而不是 x87 FPU，精度会略有不同。

相对误差非常非常小，不应破坏任何合理设计的计算。你在测试浮点相等吗？

于 2012-06-09T19:21:34.580 回答
1

我可以保证任何使用浮点数学的设备实际上都不会为数学相等的表达式提供二进制相等的结果。在您的各种平台上尝试以下 MATLAB 的等效项：

a = [repmat(1, 10000, 1); 1e16];
format long
sum(a)
sum(flipud(a))

结果：

1.000000000001000e+16
1.000000000000000e+16

加法是可交换的，所以表达式在数学上是等价的。但是顺序在浮点世界中很重要。显然，当累加值在 0 左右时，按顺序添加 10000 个 1 对于浮点应该没有问题。但是一旦浮点“浮动”到 1e16，1 太小而无法表示，所以它永远不会添加。

这是一个额外的问题：x87 FPU 在内部以扩展精度（80 位）计算。因此，如果编译器决定将中间结果保存在 FPU 中，FPU 代码会给你一个不同的答案。相反，如果它决定将中间结果溢出到堆栈，那么您将返回 64。如果它决定使用 SSE 指令进行计算，那么您将返回 64。

其他答案中建议的这些各种 MATLAB 技巧可能会让您足够接近您的问题。但是，如果您真的在对系统进行完美建模，那么您可能需要一个更可控的仿真框架。也许使用vpa, 在每一步转换为正确的位数。或者切换到 C 或 C++，特别注意优化器控制设置。

Or mathematically compute the bounds on your error based on the input scaling, and verify that your answer is always below the bound.

于 2012-06-11T13:47:28.447 回答
0

使用单精度数字，两个结果都是“相等的”（差异小于single类型的相对精度）。以下语句的计算结果为真：

max(abs( arrayfun(F,z) - F(z) )) < eps('single')

编辑

如果您真的想对此进行控制，可以尝试禁用 MATLAB 加速器以强制它对常规代码和矢量化代码使用相同的执行路径：

feature('jit', 'off')
feature('accel', 'off')
max(abs( arrayfun(F,z) - F(z) ))

feature('jit', 'on')
feature('accel', 'on')
max(abs( arrayfun(F,z) - F(z) ))

第一个/第二个的结果分别为：

ans =
     0
ans =
   5.9605e-08

显然，默认情况下加速器和即时编译器都是打开的。

于 2012-06-10T21:34:44.197 回答
0

好的。所有的答案都是正确的，但它们并不能解决问题。我想要的是数学上相等的表达式带来二进制相等的结果，无论我使用标量还是矢量算术。结果在我的 Cortex-M4、我的计算器和许多其他程序和设备中带来了 FPU，但不是 Matlab。

到目前为止，我看到的唯一解决方案是使用for-loops 而不是矢量算术。丑陋，但它的工作原理。

编辑

感谢您的回复，我们为这个问题提供了更多解决方案。

于 2012-06-11T10:04:14.763 回答