From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: LARGE roundoff error in matrix-vector multiplication?
Date: Sat, 6 Nov 2010 00:25:05 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 24
Message-ID: <ib2791$4qo$>
References: <ib1q18$51r$>
Reply-To: <HIDDEN>
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: 1289003105 4952 (6 Nov 2010 00:25:05 GMT)
NNTP-Posting-Date: Sat, 6 Nov 2010 00:25:05 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1187260
Xref: comp.soft-sys.matlab:684305

"Henry Wolkowicz" <> wrote in message <ib1q18$51r$>...
> I get the following relatively large roundoff error on matrix-vector multiplication
> >> [size(A) nnz(A) nnz(x)  normest(round(A*1e3)/1e3-A)]
> ans =
>    129   185   465   185     0
> >> norm(A(:,colPerm)*x(colPerm) - A*x)
> ans =
>   1.3523e-012
> i.e. a   A is sparse  129by185  and I am finding    A*x
> A has numbers with 2 decimal accuracy. The vector colPerm is a permutation of the integers  1:185.  By permuting the columns of A and components of x and redoing the multiplication I get an error in the 12th place. This seems like an enormous error for this type of multiplication. This does not happen with random matrices A. So perhaps there is some catastrophic 'subtraction' roundoff that happens? But getting an error in the 12th place seems enormous to me. I am trying to obtain a high accuracy algorithm and this is just checking the error at the end. So the error check has an error in it?
> Am I doing something wrong here?
- - - - - - - - - -
  You haven't told us anything about the size of values in A and x.  Your accuracy depends very much on that.  If the magnitudes of their product sums climbs up to as large as 10000 anywhere along the line during computation then one cannot expect better than rounding errors in the twelfth decimal place in your final result.  That is inherent in the concept of floating point numbers with 53-bit accuracy.  I see numbers as high as 2504 in your short example so it seems quite likely that this is the explanation for your results.

  Try this experiment.

 x = 54321*pi;
 y = x*(1/3) - x*(1/5) - x*(1/7) + x*(1/105)

Ideally y should be an exact 0 but on my computer it is off in the twelfth decimal place, and of course the explanation is that the four individual terms are mostly of the order 10^4 to 10^5 in magnitude so the inevitable round off errors are proportionally higher also.  If you left out the 54321, the error would tend to be much smaller.

Roger Stafford