罗兰谈MATLAB的艺术

将想法转化为MATLAB

一个更快更准确的总和

今天的客座博主是Christine Tobler,她是MathWorks的一名开发人员,致力于核心数字函数。
大家好!我想给你们讲一个关于舍入误差的故事,在 总和 ,以及兼容性问题。在过去的几年里,我的同事鲍比·程对 总和 让它更准确、更快,我觉得在这里讲这个故事会很有趣。

数值问题的总和

即使是一个简单的函数 总和 ,对于浮点数,求和的顺序很重要:
X = (1 + 1e-16) + 1e-16
X = 1
Y = 1 + (1e-16 + 1e-16)
Y = 1.0000
X - y
Ans = -2.2204e-16
结果是 总和 命令取决于输入相加的顺序。的 总和 MATLAB中的命令以最直接的实现开始:从第一个数字开始,一次一个地将它们相加。我在 sumOriginal 在这篇文章的底部。
但我们遇到了一个问题:对于单一精度,结果有时是非常错误的:
sumOfOnes = sumOriginal(ones(1e8, 1,“单一”))
sumOfOnes =单个16777216
这是怎么回事?问题是,对于这个数字,舍入误差大于1,我们可以通过调用 每股收益
每股收益(sumOfOnes)
Ans =单身2
因此,添加另一个1,然后四舍五入到单个精度再次返回相同的数字,因为准确的结果四舍五入到单个精度:
sumOfOnes + 1
Ans =单16777216
这一点尤其明显 总和 用于计算大型数组的平均值,因为计算的平均值可能小于输入数组的所有元素。

计算和的不同方法

现在我们已经在MATLAB中有了一个线程版本的sum,在那里我们可以并行计算一个数组的几个块的和,然后最后把它们加起来。事实证明,这个版本没有同样的问题:
numBlocks = 8;
sumInBlocks(的(1 e8, 1“单一”), numBlocks)
Ans =单个100000000
我们在MATLAB中做了这个改变 总和 即使是非线程的情况,它也解决了这种情况和其他类似的情况。但请记住,虽然这在很多情况下是一种改进,但它并不是一个完美的算法,仍然不会总是给出正确的答案。我们可以修改将输入分成的块的数量,并不是所有的结果都很好:
sumInBlocks(的(1 e8, 1“单一”), 4)
Ans =单一67108864
sumInBlocks(的(1 e8, 1“单一”), 128年)
Ans =单99999832
让我们看看另一种情况,我们不只是对数字1求和(仍然是每次都对相同的数字求和,因为这样更容易与准确的值进行比较)。我们将在这里查看结果的相对误差,这比我们得到准确的整数更相关:
X = repmat(single(3155), 54194, 1);
exactSum = 3155*54194;
numBlocks = 2.^(0:9)
numBlocks = 1×10
12 4 8 16 32 64 128 256 512
err = 0(1,长度(numBlocks));
i = 1:长度(numBlocks)
err(i) = abs(sumInBlocks(x, numBlocks(i)) - exactSum) / exactSum;
结束
重对数(numBlocks犯错)
包含(“块数”
ylabel (sumInBlocks函数中的相对错误
因此,在选择确切的块数时,肯定有一个平衡的行为,上面的图只代表一个数据集。
还有其他可能的算法来计算和,我在底部包含了一些链接。一些更复杂的问题会导致计算时间减慢,这将给没有遇到上述准确性问题的情况带来负担。
我们选择了这种“按块计算”算法,因为可以以这样一种方式进行更改,使求和变得更快。这是使用循环展开完成的:不是一次添加一个数字,而是同时将几个数字添加到单独的运行总数中。我已经包含了一个实现 sumLoopUnrolled 在下面。
顺便说一下,我并没有明确地给出选择 numBlocks 我们做了,因为我们不希望任何人依赖于它——毕竟它在未来可能会再次改变。

更改总和

在R2017b中,我们首次发布了这个新版本的 总和 ,仅适用于单精度数字的情况。对于单个数据的大型数组,所有与结果相差很大的实际问题都得到了解决,同时该函数甚至变得更快了。然而,我们也从那些对行为改变不满意的人那里得到了一些反馈;虽然新行为会带来更好的结果,但他们的代码依赖于旧行为,适应新行为是痛苦的。
虽然我们不想因为担心破坏依赖关系而陷入永远无法改进函数的困境,但我们肯定希望使更新过程尽可能轻松。一般来说,我们的目标是运行到运行的可再现性:如果一个函数用相同的输入调用两次,输出将是相同的。但是,如果机器、OS或MATLAB版本(以及其他外部因素)发生了变化,这也会在舍入误差范围内改变MATLAB返回的确切值。例如,为了性能考虑,通常会改变矩阵乘法的输出。但随着 总和 在美国,很长一段时间都没有做出任何改变,因此,越来越多的人开始依赖它的确切行为,这超出了我们的预期。
在做了这个改变之后 总和 对于单个值,我们等待几个版本来评估客户对变更的反馈。在R2020b中,我们为所有其他数据类型添加了相同的行为。这一次,我们增加了一个 发行通知 它描述了性能的改进,并提到了作为“兼容性考虑”的更改。虽然我们仍然有一些反馈,认为这是一个不方便的变化,但它比第一个情况要少。
在新的MATLAB版本中,你的一些代码是否有舍入错误的改变?你觉得发行说明中的兼容性考虑有用吗?请让我们知道 在评论中

关于准确计算和的参考文献

辅助函数

函数s = sumOriginal(x)
S = 0;
2 = 1:长度(x)
S = S + x(ii);
结束
结束
函数s = sumInBlocks(x, numBlocks)
Len = length(x)
blockLength = cell (len / numBlocks);
s = sumOriginal(x(1:blockLength));
iter = blockLength;
iter <莱恩
s = s + sumOriginal(x(iter+1:min(iter+blockLength, len));
iter = iter + blockLength;
结束
结束
函数s = sumLoopUnrolled(x)% #好< DEFNU >
numBlocks == 4的循环展开示例。为了简单起见,我们假设
% x的长度能被4整除。
注意,这种技术使用内置代码更快
%编译器标志。在这样的MATLAB代码中不一定会更快。
S1 = 0;
S2 = 0;
S3 = 0;
S4 = 0;
2 = 1:4:长度(x)
S1 = S1 + x(ii);
S2 = S2 + x(ii+1);
S3 = S3 + x(ii+2);
S4 = S4 + x(ii+3);
结束
S = s1 + s2 + s3 + s4;
结束
The MathWorks, Inc.版权所有
|

评论

如欲留言,请点击在这里登录您的MathWorks帐户或创建一个新帐户。