SuperSum, In Defense of Floating Point Arithmetic - Cleve Moler on Mathematics and Computing
 

Floating point arithmetic doesn't get the respect it deserves. Many people consider it mysterious, fuzzy, unpredictable. These misgivings often occur in discussion of vector sums. Our provocatively named [I]SuperSum[/I] is intended to calm these fears.

[B]Contents[/B]

[LIST] [*] [URL="https://www.labfans.com/bbs/#6cc58603-9b70-4a29-bc08-fae5fe3d5db0"]Ledgers[/URL] [*] [URL="https://www.labfans.com/bbs/#cc3635b6-5a76-4420-baeb-94c5b36a6573"]Checksums[/URL] [*] [URL="https://www.labfans.com/bbs/#71f16826-8ff1-4a14-bc85-146da594d56f"]Order[/URL] [*] [URL="https://www.labfans.com/bbs/#5356ff4e-b231-4d40-9404-17b0d0d882af"]Speed[/URL] [*] [URL="https://www.labfans.com/bbs/#52953716-f4f9-4698-9680-8a5d38254f1e"]Three Sums[/URL] [*] [URL="https://www.labfans.com/bbs/#f8a4c9a1-edd5-4ac3-9dca-bc7306b7de0f"]Toy Example[/URL] [*] [URL="https://www.labfans.com/bbs/#8acba015-e19b-413e-8817-1e093e31ee2d"]Second Test[/URL] [*] [URL="https://www.labfans.com/bbs/#e0e62ed1-a352-49ea-ac59-31107304d71e"]Conclusion[/URL] [/LIST]
[B]Ledgers [/B]

A [I]ledger[/I] is a list of transactions in an account. [I]Auditing[/I] the ledger involves comparing the total of the items with the change in the account balance.

If [I]v[/I] is a vector of transaction amounts, then we need to compute

sum(v) If this sum is equal to the change in the balance, then it is reasonable to assume that the ledger is correct. If not, the ledger must be examined line-by-line.

[B]Checksums [/B]

Have you ever used a checksum for a file transfer? One checksum is computed before the file is sent. After the file has been sent over a questionable channel, a second checksum is computed on the receiving end. If the two checksums agree, the transmission was probably okay. If not, the file must be sent again.

[B]Order [/B]

Floating point addition is not [I]associative[/I]. This means

(a + b) + c is not necessarily the same as

a+(b+c). So the [I]order[/I] of doing a computation is important.

Computers are deterministic devices. If the same computation is done [I]in the same order[/I] more than once, the results must be the same. If you get different sums from different runs on a fixed computer, then it must be because the order has been changed.

[B]Speed [/B]

In recent years, we have made the built-in function sum(x) faster by parallelizing it. The input vector is broken into several pieces, the sum of each piece is computed separately and simultaneously, then the partial sums are combined to give the final result. If the number and size of the pieces varies from run to run, the order varies and consequently the computed sums may vary.

[B]Three Sums [/B]

Here are three replacements for nansum(x), the version of sum(x) that skips over NaNs and infs.

[B]simplesum[/B]

Avoid the effect of blocking in the built-in sum(x).

function s = simplesum(x) % simplesum. s = simplesum(x), for vector(x). % Force sequential order for sum(x). % Skips NaNs and infs. s = 0; for k = 1:length(x) if isfinite(x(k)) s = s + x(k); end end end [B]KahanSum[/B]

This is the Kahan summation algorithm. The sum is accumulated in two words, the more significant s and the correction c. If you were able to form s + c exactly, it would be more accurate than s alone.

function s = KahanSum(x) % KahanSum. s = KahanSum(x), for vector(x). % More accurate, but slower, than sum(x). % Skips NaNs and infs. % [url]https://en.wikipedia.org/wiki/Kahan_summation_algorithm[/url]. s = 0; % sum c = 0; % correction for k = 1:length(x) if isfinite(x(k)) y = x(k) - c; t = s + y; c = (t - s) - y; s = t; end end end [B]SuperSum[/B]

I gave it a pretentious name to grab attention. Use extended precision, with enough digits to hold any MATLAB double.

function s = SuperSum(x) % SuperSum. s = SuperSum(x), for vector(x). % Symbolic Math Toolbox extended precision. % Skips NaNs and infs. % % 632 decimal digits holds every IEEE-754 double. % 632 = ceil(log10(realmax) - log10(eps*realmin)); % din = digits(632); % Remember current setting s = double(sum(sym(x(isfinite(x)),'D'))); digits(din) % Restore end [B]Toy Example [/B]

A test case here at MathWorks, known by a French-Canadian name that translates to "toy example", has a vector e2 of length 4243 and values that range from -3.3e7 to 6.8e9.

When running tests involving floating point numbers it is a good idea to set the output format to hex so we can see every last bit.

format hex load jeuTest e2 x = e2; [nansum(x) simplesum(x) KahanSum(x) SuperSum(x)] ans = 423785e8206150e2 423785e8206150e0 423785e8206150e1 423785e8206150e1 For jeuTest, Kahan summation gets the same result as SuperSum, while nansum and simplesum differ in the last bit or two.

[B]Second Test [/B]

The vector length is only three, but the third term completely cancels the first, and the second term rises from obscurity. In this situation, KahanSum is no more accurate than sum.

format hex x = [1 1e-14 -1]' [nansum(x) simplesum(x) KahanSum(x) SuperSum(x)] x = 3ff0000000000000 3d06849b86a12b9b bff0000000000000 ans = 3d06800000000000 3d06800000000000 3d06800000000000 3d06849b86a12b9b [B]Conclusion [/B]

I will leave careful timing for another day. Let's just say that in situations like jeuTest, KahanSum is probably all you need. It is usually more accurate than sum, and not much slower.

However, for complete reliability, use SuperSum. There is no substitute for the right answer.

[RIGHT][COLOR=gray][I]
[URL="javascript:grabCode_7bd084f1e3e84c86b0c419e2b2dcd8b5()"][I]Get the MATLAB code (requires JavaScript) [/I][/URL]

Published with MATLAB® R2024a
[/I][/COLOR][/RIGHT]