Accuracy of Matlab linspace and range versus Kahan summation algorithm

Question

I'm using the Matlab linspace function and the range : operator to obtain equally spaced vectors, but I'm unespectedly receiving unequally spaced numbers. My code is the following:

format long

x1 = linspace(3,5,20);
diff(x1)

x2 = 3:0.1:5;
diff(x2)

and the outputs of the vectors differences (diff) are the following:

x1    
0.105263157894737   
0.105263157894737   
0.105263157894737   
0.105263157894737
0.105263157894737
0.105263157894737
0.105263157894737   
0.105263157894737
0.105263157894737   
0.105263157894737   
0.105263157894737   
0.105263157894737
0.105263157894736   
0.105263157894737   
0.105263157894737   
0.105263157894736
0.105263157894737   
0.105263157894736   
0.105263157894737

x2
0.100000000000000   
0.100000000000000   
0.100000000000000   
0.100000000000000
0.100000000000000   
0.100000000000000   
0.100000000000000   
0.100000000000000
0.100000000000000   
0.100000000000000   
0.100000000000000   
0.100000000000001
0.100000000000000   
0.100000000000001   
0.100000000000000   
0.100000000000000
0.100000000000001   
0.100000000000000   
0.100000000000001   
0.100000000000000

To solve this problem, I'm using Kahan summation approach by the following code:

dx = 2/19;
x3 = zeros(size(x1));
x3(1) = 0;
partial_sum = 0;
c = 0.0;
for k=2:20,
    y = dx - c;
    t = partial_sum + y;
    c = (t - partial_sum) - y;
    partial_sum = t;
    x3(k) = partial_sum;
end
diff(x3)

I'm now obtaining an equispaced vector

0.105263157894737   
0.105263157894737   
0.105263157894737   
0.105263157894737
0.105263157894737   
0.105263157894737   
0.105263157894737   
0.105263157894737
0.105263157894737   
0.105263157894737   
0.105263157894737   
0.105263157894737
0.105263157894737   
0.105263157894737   
0.105263157894737   
0.105263157894737
0.105263157894737   
0.105263157894737   
0.105263157894737

However, this approach is "sequential". Is anyone aware of a vectorized or parallel implementation of the Kahan summation approach to improve the efficiency or porting it to a (CUDA) parallel reduction?

Thank you in advance.

worth a look: [How does the COLON operator work?](http://web.archive.org/web/20120213165003/http://www.mathworks.com/support/solutions/en/data/1-4FLI96/index.html?) (retrieved using Wayback Machine, MathWorks removed page for some reason!) — Amro, Apr 19 '13 at 09:03

score 4 · Answer 1 · answered Apr 19 '13 at 09:03

4

I think that you are misunderstanding

(a) the use of linspace. Your call

linspace(3,5,20)

returns 20 equally spaced numbers starting at 3 and ending at 5. Try

linspace(3,5,21)

instead. Note that the lengths of your vectors x and x2 differ by one.

and (b)

the subtleties of floating-point arithmetic. This topic is covered extensively on SO and I won't reproduce any of the many good explanations already to be found.

I think that what you are seeing are numbers which are equally spaced enough for most numerical applications.

answered Apr 19 '13 at 09:03

High Performance Mark

77,191
7
105
161

you are right, but I think this was an off-by-one error. Plus the example used was not the best to show the difference. Try this instead: `(0:0.1:1) - linspace(0,1,11)` – Amro Apr 19 '13 at 09:06
I think that we are all overlooking the main question: `Kahan summation`. – fpe Apr 19 '13 at 09:11
Well, possibly, and my 'answer' is, I acknowledge, more of an extended comment than an answer to the (part of the) question about Kahan summation. But, if what OP means is 'I haven't quite figured out linspace and : so I'm going to implement Kahan summation' my answer is 'figure out linspace and : first'. But hey ho, whatever dude. – High Performance Mark Apr 19 '13 at 09:13
@Amro Why should I try `(0:0.1:1) - linspace(0,1,11)`? – Vitality Apr 19 '13 at 12:36
@fpe Yes. More precisely: vectorization/parallelization of Kahan summation method or similar methods. – Vitality Apr 19 '13 at 12:36
@JackOLantern: it was just to show that `linspace` and `colon` do not always produce the same thing. Not exactly what you were asking I guess.. – Amro Apr 19 '13 at 13:11

fpe · Answer 2 · 2013-04-19T09:13:37.397

I guess that this link is gonna be helpful:

http://www.mathworks.com/matlabcentral/fileexchange/26800

You will find the XSum function, where you have the chance to switch and test different summation methods:

Double: A thread-safe implementation of Matlab's SUM. At least in Matlab 2008a to 2009b the results of the multi-threaded SUM can differ slightly from call to call. 50% faster than single-threaded SUM (MSVC++ 2008), equivalent accuracy.
Long: Accumulated in a 80 bit long double, if the compiler support this (e.g. LCC v3.8). 3.5 more valid digits, 40% slower.
Kahan: The local error is subtracted from the next element. 1 to 3 more valid digits, 10% slower.
Knuth: As if the sum is accumulated in a 128 bit float: about 15 more valid digits. About the same speed (MSVC++2008 compiler). This is suitable for the most real world problems.
Knuth2: 30 more valid digits as if it is accumulated in a 196 bit float. 60% slower.
KnuthLong: As Knuth, but using long doubles to get about 21 more valid digits, if supported by the compiler. 2.5 times slower.

Thank you very much for your answer. Could you please shortly comment on vectorization/parallelization? You talk about multi-threaded sum, at least for the "Double" method. Is that already "multi-thread-ready"? What about the other methods? Thanks. — Vitality, Apr 19 '13 at 12:29

Accuracy of Matlab linspace and range versus Kahan summation algorithm

2 Answers2