# bug #54567: median(uint8) gives bad results

 Submitted by: None Submitted on: Fri 24 Aug 2018 03:50:40 PM UTC Category: Octave Function Severity: 3 - Normal Priority: 5 - Normal Item Group: Incorrect Result Status: Patch Submitted Assigned to: None Originator Name: Etienne de Foras Originator Email: -email is unavailable- Open/Closed: Open Release: 4.4.1 Operating System: Any

You are not logged in

## Discussion

Mon 27 Aug 2018 11:44:42 PM UTC, comment #16:

I implemented many of the ideas in the various comments in to a patch which is attached.  The case for int64 is not ideal since it warns when there is potential loss of accuracy, but should just calculate an accurate result.  The code is tricky enough that maybe Dan's suggestion of using mean() would be better.  First, fix mean() to handle every single case and then let median call mean.

One possibility is to use the algorithm that Dan suggested in comment #11.  One could use permute ahead of the algorithm in order to guarantee that the correct dimension was lined up the way the algorithm expects.  After calculating the median one would then use ipermute to reverse the change.

Rik <rik5>
Mon 27 Aug 2018 12:19:28 AM UTC, comment #15:

Dan: Right. Just thought that this might be of interest,
and, of course, it is a "signal-processing" technique.

Michael Godfrey <godfrey>
Sun 26 Aug 2018 09:21:07 PM UTC, comment #14:

The Mohanty algorithm is actually for "running" median, or another common terminology, "moving median".  It keeps a sorted linked list for which it can easily drop a sample when the window moves past that sample, and add a sorted sample in the list when that sample moves into the window.  This would be something good for the signal-processing toolbox if a moving median doesn't exist already.  (Of course, Octave scripting of linked lists isn't very efficient, so written in C++ is preferred.)

Octave's current median() with sorting is pretty efficient.  The mean() function definitely needs some work, though.  There probably is not a lot of coding, but because of its importance the routine needs attention to detail which could take time.

Dan Sebald <sebald>
Sun 26 Aug 2018 01:08:39 PM UTC, comment #13:

I added comment #6 just for information. A median algorithm for
very large datasets (as, for example, the LIGO data)
is fully discussed in:

This reference includes the C-code for the algorithm.

Also, the current Octave mean() just uses sum()/n which is
not robust nor efficient. Improving this might be a good
project.

Michael Godfrey <godfrey>
Sun 26 Aug 2018 07:29:38 AM UTC, comment #12:

BTW, we should at BIST tests for all these cases.

Also, we know mean() probably has the issue we are trying to avoid:

octave:253> mean([intmin('int64')+5 intmax('int64')-5],'native')
ans = 0
octave:254> mean([int64(-1) int64(0)],'native')
ans = -1

I.e., the large magnitude int64 numbers have a loss of precision in mean() so the result doesn't match the formula of Comment #11.

Dan Sebald <sebald>
Sun 26 Aug 2018 07:25:46 AM UTC, comment #11:

KT, the following seems to work for int64:

octave:234> dim = 2;
octave:235> xs64 = [1 6 11 16; -20 -5 14 19; intmin('int64') intmin('int64')+5 intmax('int64')-5 intmax('int64')]
xs64 =

1                     6                    11                    16
-20                    -5                    14                    19
-9223372036854775808  -9223372036854775803   9223372036854775802   9223372036854775807

octave:236> f = (sum (sign(x), dim) != 0)
f =

1
0
0

octave:237> f .* x(:,1) + (x(:,2) + !f .* x(:,1) - f .* x(:,1)) / 2
ans =

9
5
-1

Notice there is one small issue with the fact I had to use the colon indexing, meaning that I'm assuming dim==2.  That means you'll have to have an if-else statement based on "dim" with slightly different formulas for dim=1 versus dim=2.  That's unfortunate.  I wonder if we should devise an indexing type of function that essentially extracts a sub-matrix from a matrix the way so many functions can act on "dim", e.g., nth_element(A,k,dim), sum(A,dim), mean(dim).

Also, if this approach proves more efficient than all the min/maxing of the uint64 formula in Comment #10, perhaps use of logical arrays could be used for unit64 as well.

Dan Sebald <sebald>
Sun 26 Aug 2018 05:22:16 AM UTC, comment #10:

KT, here's a test of the unsigned 64-bit integer using min/max functions:

octave:110> x = nth_element (x64, 2:3, 2);
octave:111> min (x, [], 2) + (max (x, [], 2) - min (x, [], 2)) / 2
ans =

9
12
18446744073709551608

[The following does not work generally]
octave:112> x(:,1) + (x(:,2) - x(:,1)) / 2
ans =

9
12
18446744073709551608

octave:113> x(:,2) + (x(:,1) - x(:,2)) / 2
ans =

11
14
18446744073709551610

I tested what happens when the unsigned "wrap" occurs.  Well, Octave is programmed to not wrap.  I was curious what happens with wrapping because sometimes one can devise a formula that allows integer wrapping so long as it is guaranteed to unwrap.  Since Octave truncates rather than wraps, there's no tricks that can be done.

OK, so let's think about an int64 formula then...

Dan Sebald <sebald>
Sun 26 Aug 2018 04:45:47 AM UTC, comment #9:

That's fine with me.

On the matter of two numbers, for what it is worth, it seems the floating point conversion (then sum, then divide) is the fastest.  Consider

octave:34> x = uint8(255*ones(100000,1));
octave:35> y = uint8(254*ones(100000,1));
octave:36> tic; zfloat = mean([x y], 2, "native"); toc
Elapsed time is 0.00436807 seconds.
octave:37> tic; zintwrong = bitshift(x + y + 1, -1); toc
Elapsed time is 0.00426006 seconds.
octave:38> tic; zintright = bitshift(uint32(x) + uint32(y) + 1, -1); toc
Elapsed time is 0.00503707 seconds.

So the additional conversion to uint32, rounding and division-by-bitshift comes out slower than converting to float and doing the math.  Even when x and y vectors are not converted, it is still on the order of converting to float first.  I would think that in assembly code that the unsigned int conversion might be "free", i.e., just a case of a different machine instruction that treats the higher register bits as an extension from 8 bits to 32 bits.  I suppose what the floating point instructions have going for them is inherent rounding.

Dan Sebald <sebald>
Sat 25 Aug 2018 06:03:16 PM UTC, comment #8:

Even though it makes for clear code to call mean(), we should not do it.  This is the mean of only two numbers and the overhead involved in calling an Octave function is significant.  For core Octave we do prioritize speed and efficiency to some degree over readability of code.  So the calculation should be be done in median.m, and probably needs an exception case for (u)int64 values.  I can review any amended changeset.

Rik <rik5>
Sat 25 Aug 2018 04:30:25 PM UTC, comment #7:

Nice article on mean computation...  This is why I wondered if it were better to put this case of 64-bit average in the mean() script and call mean() inside median() instead.  Could there be an internal library routine that is stable that making a built-in mean() makes sense?

Still, here for median() it is only two numbers we are averaging, plus the formula in the reference pertains to float values, not integer "wrap" effects.  Consider that the formula from the reference becomes

mu_2 = mu_1 + 1/2 (x_2 - mu_1)
= x_1  + 1/2 (x_2 - x_1)

which is the same formula as in Comment #4.

KT, keep thinking about the formulas; it shouldn't be too difficult to implement conditionals with scripts.  They slow things down, but I'd prefer being correct than have a warning message.  For example, the formula of Comment #4 could be implemented as

retval = min (x, dim) + (max (x, dim) - min (x, dim))/2;

Maybe there is a way to generalize the formula from the reference for integers.  The issue with integers is the loss of the fractional portion when doing the intermediate mean in the recursion.  Perhaps there is a way to accumulate those lost bits as some sort of residual and then at the very end of the computation adjust the final computation with the residual divided by N.

Dan Sebald <sebald>
Sat 25 Aug 2018 01:53:59 PM UTC, comment #6:

Numerically stable mean value algorithms have been around
for a long time. The one shown in the URL below dates from
the 1960's:

The original references should be easy to find.

Michael Godfrey <godfrey>
Sat 25 Aug 2018 12:12:20 PM UTC, comment #5:

Good catch with the (u)int64 data type.

I think the second version of comment #2 will deal with averaging two numbers for the even case with (u)int8-32, single, double, and logical.  For (u)int64, we might have to introduce an exception

## Perform the averaging in double precision to avoid intermediate
## overflows of finite data types like "int8".  For "int64" and
## "uint64" double precision is insufficient and we have to
## perform the division by two beforehand.

x = nth_element (x, k:k+1, dim);
if ((isa (x, "int64") || isa (x, "uint64"))
&& (any (x >= flintmax () / 2)))
warning ("median: inaccurate result due to large integer values.")
retval = sum (x / 2, dim);
else
retval = cast (sum (x, dim) / 2, class (x));
end

The formula of comment #3 and comment #4 is hard to implement regarding multi-dimensional arrays.  If you found a more elegant way, I would be very happy ;-)

Kai Torben Ohlhus <siko1056>
Sat 25 Aug 2018 12:53:29 AM UTC, comment #4:

"For uint64, it's just meanXY = X + (Y-X)/2."

Actually, need to ensure that we pick Y >= X or there is wrap-around that way.

Dan Sebald <sebald>
Fri 24 Aug 2018 11:51:03 PM UTC, comment #3:

This is one of those I'd classify as "think about for a while".  I did notice the "logical" case, both within median() and mean().  However, I din't think why it is there, just that mean() with the follow-up check within median() should still work fine.  Of course, in the alternative changeset there is an extra call as well, to cast().

Yes, mean() is an extra function call of overhead, as mean() in turn calls sum().  For large data this overhead should be diminished.  There's nothing recursive or loop-like about these routines, so the greater the input matrix size, the more CPU goes toward the inner matrix operations written in C++.

I'm fine with dropping the "native" and then use cast().  I was thinking more along the lines of there being just one location where the sophistication of mean() is implemented.  The one thing to be careful of, though, is something we've run into with variable indexing, i.e., a double has 52 bit mantissa or significand.  Hence, there is this 12-bit range where int-64 numbers lose resolution when computed with double values.

The general case of mean() means this is complicated.  We are summing together unknown N numbers, so we need ceil(log2(N)) bits overhead to ensure there is no arithmetic overflow.  But in the case of median() it is just two numbers that involve mean.  So, should we consider the special cases of median of int64 and uint64?

For int64, given two numbers X and Y, if X and Y are of different sign, then meanXY = (X+Y)/2.  If X and Y are the same sign, then meanXY = (2X - X + Y)/2 = X + (Y-X)/2, or something close ensuring rounding is appropriate.

For uint64, it's just meanXY = X + (Y-X)/2.

Keep in mind that this needs to be for vector/matrix data, so the test above would involve logical arrays in some way.

Dan Sebald <sebald>
Fri 24 Aug 2018 10:54:14 PM UTC, comment #2:

Matlab R2018a returns "140" in both cases of the original comment.

The patch file #44857 of comment #1 works

-    retval = sum (nth_element (x, k:k+1, dim), dim, "native") / 2;
+    retval = mean (nth_element (x, k:k+1, dim), dim, "native");

But I think it makes things more complicated and maybe slower for large data?!  Basically, mean() would perform the computation with double precision and finally convert the result to the input data type.  Why not doing the very same without additional function calls in one line and tidy up the nasty "logical" exception:

-    retval = sum (nth_element (x, k:k+1, dim), dim, "native") / 2;
-    if (islogical (x))
-      retval = logical (retval);
-    endif
+    ## Perform the averaging in double precision to avoid intermediate
+    ## overflows of finite data types like "int8".
+    retval = cast (sum (nth_element (x, k:k+1, dim), dim) / 2, class (x));

I thought about other solutions from the "real plain dump" integer world (like C/C++) but these would not apply for Octave... This I am afraid to stay to this solution.

Kai Torben Ohlhus <siko1056>
Fri 24 Aug 2018 07:22:27 PM UTC, comment #1:

The question is whether Matlab does the arithmetic using the intrinsic variable type.  The median() routine is supposed to return the same variable type, but that doesn't necessarily mean the average of the center values should be computed using that type.  Why should it, I guess?  In fact, mean() appears to do its arithmetic not using native, but instead converts it afterward:

octave:12> u = [1 140 140 255];
octave:13> u8 = uint8(u);
octave:14> mean(u)
ans =  134
octave:15> mean(u8)
ans =  134

So, attached a patch to change the behavior for median.  Basically it uses mean() instead of sum() with the assumption that the proper averaging behavior is already baked into mean().  I added a couple more tests as well.

Dan Sebald <sebald>
Fri 24 Aug 2018 03:50:40 PM UTC, original submission:

If I define:

u= [1 140 140 255];
u8=uint8(u);

then median(u)=140 (correct)
but  median(u8)=128 (incorrect, should be 140)

The bug disappear with odd size data, the median with even size data is computed as (140+140)/2, but with saturation to 255 due to uint8 type and then divided by two ->127 (or 128 with rounding).
We may use the type int32 or double for the sum and division by 2 and then convert back to uint8

(Thanks for GNU Octave !)

Anonymous

## Attached Files

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:

Comment:

Attached Files
file #44882:  54567.cset added by rik5 (2KiB - application/octet-stream)

## Dependencies

Depends on the following items: None found

Items that depend on this one: None found

Carbon-Copy List
• -email is unavailable- added by rik5 (Posted a comment)
• -email is unavailable- added by godfrey (Posted a comment)
• -email is unavailable- added by siko1056 (Posted a comment)
• -email is unavailable- added by sebald (Updated the item)
• -email is unavailable- added by None (Submitted the item)
•

Do you think this task is very important?
This task has 0 encouragements so far.

Only project members can vote.