Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
accumulation error caused by diag and sum

Subject: accumulation error caused by diag and sum

From: lily

Date: 6 Jun, 2012 09:07:08

Message: 1 of 5

I found the following situation, combined "diag" and "sum" commands will cause accumulation error. Then how to eliminate this condition? Thanks.

B=1e+10*rand(25);
B=B-diag(diag(B));
sumB=sum(B);
C=B;
C=C+diag(-sum(B));
sumC=sum(C);

err=sumC(abs(sumC)>1e-5); %what I need is the matrix C, sum(C) should be all zeros.
tf1=any(diag(B)); %tf1=0 means code B=B-diag(diag(B)) produces no error.
tf2=any(sumB~=sum(B)); %tf2=0 means "sum" produces no error.

>>err = % how to eliminate accumulation error?
  1.0e-004 *
    0.1144 -0.1907
>>tf1 = 0
>>tf2 = 0

Subject: accumulation error caused by diag and sum

From: John D'Errico

Date: 6 Jun, 2012 13:52:07

Message: 2 of 5

"lily" wrote in message <jqn6jr$gdp$1@newscl01ah.mathworks.com>...
> I found the following situation, combined "diag" and "sum" commands will cause accumulation error. Then how to eliminate this condition? Thanks.
>
> B=1e+10*rand(25);
> B=B-diag(diag(B));
> sumB=sum(B);
> C=B;
> C=C+diag(-sum(B));
> sumC=sum(C);
>
> err=sumC(abs(sumC)>1e-5); %what I need is the matrix C, sum(C) should be all zeros.
> tf1=any(diag(B)); %tf1=0 means code B=B-diag(diag(B)) produces no error.
> tf2=any(sumB~=sum(B)); %tf2=0 means "sum" produces no error.
>
> >>err = % how to eliminate accumulation error?
> 1.0e-004 *
> 0.1144 -0.1907
> >>tf1 = 0
> >>tf2 = 0

Welcome to the world of floating point arithmetic.

You can't eliminate all such errors. Not possible here.

In fact, the order in which you perform even simple
addition and subtraction operations will affect your
result when done in floating point arithmetic.

Instead, learn to work with it. Learn to make your
computations insensitive to such problems. Classes
and books on numerical methods teach you just
that.

John

Subject: accumulation error caused by diag and sum

From: Jan Simon

Date: 6 Jun, 2012 15:22:09

Message: 3 of 5

Dear lily,

You can use a sum with error compenstaion:
  http://www.mathworks.com/matlabcentral/fileexchange/26800-xsum
Here the intermediate results are stored with a higher precision. This catchs problems in e.g.:

  x = [1e17, 1, -1e17]
  sum(x)
  % >> 0
  XSum(x)
  % >> 1

XSum can use 128 or 196 bits for the accumulation. This helps, but it does not solve the problem of the limited precision:

  XSum([1e100, 1e-100, -1e100])

This will still be not accurate.

Kind regards, Jan

Subject: accumulation error caused by diag and sum

From: lily

Date: 7 Jun, 2012 05:51:07

Message: 4 of 5

Thanks for John and Jan's advice, it helps me a lot. Now the code is edited as following:

B=1e+10*rand(20);
B=B-diag(diag(B));
sumB=sum(B);
C=B;
C=C+diag(-sum(B));
% improve error tolerance and compensate error
% outside of the program,abs(sum(C))<1e-4 is considered as sum(C)'s all elements be zeros.
sumC=sum(C);
ids=find(abs(sumC)>=1e-5);
errs=sumC(ids);
for i=1:length(ids)
    id=ids(i);
    C(id,id)=C(id,id)-errs(i);
end
sumC=sum(C);
ids=find(abs(sumC)>=1e-5);
errs2=sumC(ids);
>>errs =
  1.0e-004 *
  Columns 1 through 9
    0.1717 -0.1431 -0.2098 0.3052 0.1574 0.1907 0.1764 0.1335 -0.1097
  Column 10
    0.1144
>>errs2 =
 -1.2875e-005

Subject: accumulation error caused by diag and sum

From: Bruno Luong

Date: 7 Jun, 2012 06:53:07

Message: 5 of 5

"lily" wrote in message <jqn6jr$gdp$1@newscl01ah.mathworks.com>...
> I found the following situation, combined "diag" and "sum" commands will cause accumulation error. Then how to eliminate this condition? Thanks.
>
> B=1e+10*rand(25);
> B=B-diag(diag(B));
> sumB=sum(B);
> C=B;
> C=C+diag(-sum(B));
> sumC=sum(C);
>
> err=sumC(abs(sumC)>1e-5); %what I need is the matrix C, sum(C) should be all zeros.

You can't make a reliable code based on SUM. SUM order and internal precision is not specified by TMW. So the results might vary from computer-to-computer, MATLAB versions.

To make a good numerical code, you have to chose an algorithm that is robust on lost of precision due to floating point calculation.

Bruno

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us