Which one is the right answer, sum(w(:)) or sum(sum(w))?

Hi, all

Just as title. When process certain matrix, like the one in attachment, I found different method returns different result, which one should be used?

>> sum(w(:)) - sum(sum(w))

ans =

single

-3.8147e-06

Since I am quite new to MATLAB, can anyone tell me, is there any other function or Best Practical method to solve such problem?

I know, the real cause of this problem is due to the limitation of bit length on hardware, but if the result is highly depend on the order of matrix elements, I have to wondering neither one of the methods can get the right answer.

This issue is related to sum(), but also influence related functions like mean(), etc. Since all these are daily used element blocks, I thinks, it worth to be clearified in more detailed.

After all, accuracy is critical to mathmatic, and the truth is only one. The aim of MATLAB is to make user focus on the beauty of math but not the detail of hardware, isn't it?

Waiting for any comments!

Regards,

Liu Lan

Answer 1

[Walter Roberson]

For sum(sum) and sum((:)) and sum('all'), provided that the array is large enough, one or more calls will be made to a high speed mathematics library that is automatically multi-threaded. The data is divided into sections and each section is totaled and then the partial sums are added. When you use sum(sum) versus sum((:)) the boundaries of the partial sums may be different, and due to floating point roundoff, you can get different results.

Mathworks does not provide a "compensated" additional function designed to remove the roundoff error.

Answer 2

[Matt J]

sum(sum(w)) will generally be slower, because it requires 2 function calls. Also, it only works for 2D arrays whereas sum(w(:)) or sum(w,'all') work for w of any dimension.

There is no way you can optimize the accumulation of floating point errors through any one of the 3 calling strategies. The error is influenced by the way the summation is split into parallel jobs and distributed to your multicore processor. This is done in ways that are hidden from you as the user, and which will be hardware and OS dependent.

[Walter Roberson]

w(:) turns out to be a function call for matlab purposes. If w is marked as complex but all of the complex components are zero or negative zero, then w(:) will drop the complex component. The data pointer for (:) will be different than the original one for complex values.

reshape(w, [], 1) on the other hand is treated specially and will retain any all-zero complex component. The data pointer will not be changed: reshape creates a new dimension descriptor pointing to the same data.

Answer 3

[Bruno Luong]

"Since I am quite new to MATLAB, can anyone tell me, is there any other function or Best Practical method to solve such problem?"

The best pratice IMO is using

sum(w,'all')

if your MATLAB can support it. I believe.

It seems the result is the same than sum(w(:)) but avoid the tiny slowness of (:) for complex data. NOTE that Walter warns about that in his comment, however I'm not 100% sure the Excution engine is smart enough to bypass explicit call of (:) in the statement sum(w(:)), but let us assume it is not.

The order has been studied in depth by TMW and they come up with relative robust but fast way to implement sum for generic data.

The most accurate otder depends on your data (oscillante data like Bessel function in long rage is horrible to compute caccuretely), and the good order is very difficult to know and very expansive to estimate. I don't believe sum is designed for that purpose, and there is few implementations on FEX when I look a while ago.

But your ^rincipal duty is design the overall numerical method that should be less sensitive to truncation error of sum, mean etc... generally anything related to limited floating point precision. If not your algorithm is a bad one.

[Lan]

Does FEX means 3rd parties extension package?

[Bruno Luong]

FEX is File exchange. I retrieve a package call AccSum in 2010, based on the enclosed articles.

As it is not my code I cannot post it here.

[Lan]

Many thanks for your input!

Answer 4

[Lan]

Hi,

First of all, I want to say thanks to all the people response to my question! And I am quite impressived by the warmly activities on this forum!

After that, I have to say, none of you have caught the key to my question:

Is there any method to make sum() stick to the Right Answer of a 'single' type matrix?

The definition of Right Answer here should be the unique 'single' value, rounded from Exactly Arithmatic summation of all members in specified matrix.

Anyway, from the discussion, I think I might get the answer: there is no such build-in function in MATLAB can achieve that purpose, at least on 'single' type, isn't it?

I don't think this problem is not solvable, even using parallel computing on multicore. Because, I have try the same script on different type of CPUs, i5 dual-core, i7 quad-core, Xeon, including online version. The result error is exactly same. Which means, the compute logic is still under control, and the result is by "Design".

From my understanding, make function run faster is good, but usually, it is depend on the accuracy as premise. Take a step back, why not provide a build-in "slow but precise" version of sum() just with the proper counting order?

Some guy said it is not necessory to keep such precise. I have to say, it only depends on the application. If the error is introduced in early stage of the calculation, it can be easily expanded to 0.1% different in result.

By the way, I think it's a common way to use mean() to remove bias of a matrix in the first step, isn't it? This will introduce such kind of error.

Others might say, if precise is important, why not using 'double'. Yes, for my case, I can use this method. But does it means the 'single' type is useless in computing, except for more compact store?

Meanwhile, I just wondering, with following information, the same problem might happened, even on 'double' case.

https://www.mathworks.com/help/matlab/release-notes.html?category=performance&rntext=&startrelease=R2020b&endrelease=R2020b&groupby=release&sortby=descending

I really hope this kind of issue can be handled more elegent by Mathworks.

Regards,

Liu Lan

[Walter Roberson]

Is there any method to make sum() stick to the Right Answer of a 'single' type matrix?

We already told you NO.

@Jan has a File Exchange contribution that implements several different extended precision summation methods; see https://www.mathworks.com/matlabcentral/fileexchange/26800-xsum

Any of the methods in his description that talk about being "equivalent" to a certain number of decimal places, will not give results "rounded from Exactly Arithmatic summation"

Several years ago, I read a paper in which some numeric analysis people demonstrated a family of functions such that for any given binary floating point precision (e.g., could be 10000 binary digits), that a function can be constructed such that the evaluation of the function will, in an infinite number of places, round incorrectly compared to mathematical analysis.

Is it theoretically possible to get correct summation if you stick to the input values being IEEE 754 Double Precision ? Maybe; it might require aligning everything in a 2098 (or possibly more) bit number. (abs(log2(eps(realmin))) + log2(realmax) = 2098 - you might be adding the largest possible number and the smallest possible number, and you cannot lose any bits because the expression might be

SUM([realmax(), eps(realmin), -realmax())

and the result needs to be eps(realmin) exactly.

One thing you can do is use

sum(sym(YourDoubleMatrix, 'f'))

as the sym() part will convert each value to its internal rational value, such as 884279719003555/281474976710656 being the internal representation for the closest representable value to π . The symbolic toolbox can then add the rationals, and you can do appropriate rounding afterwards.

[Bruno Luong]

@Lan "The definition of Right Answer here should be the unique 'single' value, rounded from Exactly Arithmatic summation of all members in specified matrix.

Anyway, from the discussion, I think I might get the answer: there is no such build-in function in MATLAB can achieve that purpose, at least on 'single' type, isn't it?"

There is NOT such floating point algorithm based on native sum of 2 floating numbers, MATLAB or not.

[Walter Roberson]

https://dl.acm.org/doi/10.1145/3544488 seems to imply that there are known algorithms

[Lan]

I am really sorry for my dull of can not understanding the implict expression in the first time. And many thanks for your time for make it clearify again!