http://www.mathworks.com/matlabcentral/newsreader/view_thread/295702
MATLAB Central Newsreader  LARGE roundoff error in matrixvector multiplication?
Feed for thread: LARGE roundoff error in matrixvector multiplication?
enus
©19942014 by MathWorks, Inc.
webmaster@mathworks.com
MATLAB Central Newsreader
http://blogs.law.harvard.edu/tech/rss
60
MathWorks
http://www.mathworks.com/images/membrane_icon.gif

Fri, 05 Nov 2010 20:39:04 +0000
LARGE roundoff error in matrixvector multiplication?
http://www.mathworks.com/matlabcentral/newsreader/view_thread/295702#793596
Henry Wolkowicz
I get the following relatively large roundoff error on matrixvector multiplication<br>
<br>
>> [size(A) nnz(A) nnz(x) normest(round(A*1e3)/1e3A)]<br>
ans =<br>
129 185 465 185 0<br>
>> norm(A(:,colPerm)*x(colPerm)  A*x)<br>
ans =<br>
1.3523e012<br>
<br>
i.e. a A is sparse 129by185 and I am finding A*x<br>
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?<br>
Am I doing something wrong here?

Fri, 05 Nov 2010 22:16:03 +0000
Re: LARGE roundoff error in matrixvector multiplication?
http://www.mathworks.com/matlabcentral/newsreader/view_thread/295702#793611
Bruno Luong
"Henry Wolkowicz" <hwolkowicz@uwaterloo.ca> wrote in message <ib1q18$51r$1@fred.mathworks.com>...<br>
> I get the following relatively large roundoff error on matrixvector multiplication<br>
> <br>
> >> [size(A) nnz(A) nnz(x) normest(round(A*1e3)/1e3A)]<br>
> ans =<br>
> 129 185 465 185 0<br>
> >> norm(A(:,colPerm)*x(colPerm)  A*x)<br>
> ans =<br>
> 1.3523e012<br>
> <br>
> i.e. a A is sparse 129by185 and I am finding A*x<br>
> 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?<br>
> Am I doing something wrong here?<br>
<br>
<a href="http://en.wikipedia.org/wiki/Condition_number">http://en.wikipedia.org/wiki/Condition_number</a><br>
<br>
Your matrix is badly scaled.<br>
<br>
Bruno

Fri, 05 Nov 2010 23:04:04 +0000
Re: LARGE roundoff error in matrixvector multiplication?
http://www.mathworks.com/matlabcentral/newsreader/view_thread/295702#793622
Henry Wolkowicz
Bruno ... your reply does not make sense to me. I am doing matrixvector multiplication. What does an illconditioned matrix (if it is or not) have anything to do with matrix multiplication? <br>
In fact, this can be considered to be the inner product of two vectors (looking at just one row of the matrix vector multiplication) with only 6 numbers.<br>
<br>
1.0e+003 *<br>
<br>
0.001000000000000 0.000000000000000 0.001000000000000 0.000000000000000<br>
0.000400000000000 0.158000000000000 0.001700000000000 0.158000000000000<br>
0.000400000000000 0.158000000000000 0.001500000000000 2.504373333333334<br>
0.001500000000000 2.504373333333334 0.000400000000000 0.158000000000000<br>
0.001700000000000 0.158000000000000 0.000400000000000 0.158000000000000<br>
0.001700000000000 0.158000000000000 0.001700000000000 0.158000000000000<br>
<br>
when I multiply (inner product) the first two vectors and the last two I get<br>
<br>
valrow3 =<br>
<br>
3.092960000000001e+003<br>
<br>
<br>
valrow3perm =<br>
<br>
3.092960000000001e+003<br>
<br>
<br>
ans =<br>
<br>
4.547473508864641e013<br>
<br>
<br>
where the last 'ans' is the subtraction of the two values of the inner products. The numbers look exactly the same to me ... but the subtraction is still of the order<br>
of 1e12<br>
<br>
<br>
<br>
<br>
"Henry Wolkowicz" <hwolkowicz@uwaterloo.ca> wrote in message <ib1q18$51r$1@fred.mathworks.com>...<br>
> I get the following relatively large roundoff error on matrixvector multiplication<br>
> <br>
> >> [size(A) nnz(A) nnz(x) normest(round(A*1e3)/1e3A)]<br>
> ans =<br>
> 129 185 465 185 0<br>
> >> norm(A(:,colPerm)*x(colPerm)  A*x)<br>
> ans =<br>
> 1.3523e012<br>
> <br>
> i.e. a A is sparse 129by185 and I am finding A*x<br>
> 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?<br>
> Am I doing something wrong here?

Fri, 05 Nov 2010 23:26:04 +0000
Re: LARGE roundoff error in matrixvector multiplication?
http://www.mathworks.com/matlabcentral/newsreader/view_thread/295702#793626
Bruno Luong
"Henry Wolkowicz" <hwolkowicz@uwaterloo.ca> wrote in message <ib22h4$9bp$1@fred.mathworks.com>...<br>
> Bruno ... your reply does not make sense to me. I am doing matrixvector multiplication. What does an illconditioned matrix (if it is or not) have anything to do with matrix multiplication? <br>
<br>
You are right, sorry somehow I though you solve linear system.<br>
<br>
<br>
> <br>
> valrow3 =<br>
> <br>
> 3.092960000000001e+003<br>
> <br>
> <br>
> valrow3perm =<br>
> <br>
> 3.092960000000001e+003<br>
> <br>
> <br>
> ans =<br>
> <br>
> 4.547473508864641e013<br>
> <br>
> <br>
> where the last 'ans' is the subtraction of the two values of the inner products. The numbers look exactly the same to me ... but the subtraction is still of the order<br>
> of 1e12<br>
> <br>
<br>
4.547473508864641e013 is eps(3.092960000000001e+003), so what you get is simply at a limit of floatingpoint accuracy. Is that surprise you? Not me.<br>
<br>
Bruno

Fri, 05 Nov 2010 23:48:05 +0000
Re: LARGE roundoff error in matrixvector multiplication?
http://www.mathworks.com/matlabcentral/newsreader/view_thread/295702#793630
James Tursa
"Henry Wolkowicz" <hwolkowicz@uwaterloo.ca> wrote in message <ib22h4$9bp$1@fred.mathworks.com>...<br>
> Bruno ... your reply does not make sense to me. I am doing matrixvector multiplication. What does an illconditioned matrix (if it is or not) have anything to do with matrix multiplication? <br>
> In fact, this can be considered to be the inner product of two vectors (looking at just one row of the matrix vector multiplication) with only 6 numbers.<br>
> <br>
> 1.0e+003 *<br>
> <br>
> 0.001000000000000 0.000000000000000 0.001000000000000 0.000000000000000<br>
> 0.000400000000000 0.158000000000000 0.001700000000000 0.158000000000000<br>
> 0.000400000000000 0.158000000000000 0.001500000000000 2.504373333333334<br>
> 0.001500000000000 2.504373333333334 0.000400000000000 0.158000000000000<br>
> 0.001700000000000 0.158000000000000 0.000400000000000 0.158000000000000<br>
> 0.001700000000000 0.158000000000000 0.001700000000000 0.158000000000000<br>
> <br>
> when I multiply (inner product) the first two vectors and the last two I get<br>
> <br>
> valrow3 =<br>
> <br>
> 3.092960000000001e+003<br>
> <br>
> <br>
> valrow3perm =<br>
> <br>
> 3.092960000000001e+003<br>
> <br>
> <br>
> ans =<br>
> <br>
> 4.547473508864641e013<br>
> <br>
> <br>
> where the last 'ans' is the subtraction of the two values of the inner products. The numbers look exactly the same to me ... but the subtraction is still of the order<br>
> of 1e12<br>
<br>
You can try these different methods to see numerical differences just by reordering the calculations:<br>
<br>
Pick of two of your vectors for the dot product calculation, call them A and B, and then try these methods:<br>
<br>
dot(A,B)<br>
A.'*B<br>
mtimesx(A,'t',B,'blas')<br>
mtimesx(A,'t',B,'loops')<br>
mtimesx(A,'t',B,'loopsomp')<br>
<br>
The last one requires an OpenMP capable compiler and a large size for A and B to show any difference from the 2nd to last one. MTIMESX can be found here:<br>
<br>
<a href="http://www.mathworks.com/matlabcentral/fileexchange/25977mtimesxfastmatrixmultiplywithmultidimensionalsupport">http://www.mathworks.com/matlabcentral/fileexchange/25977mtimesxfastmatrixmultiplywithmultidimensionalsupport</a><br>
<br>
Based on some of my testing, dot(A,B) is the slowest and least accurate. The 'loops' or 'loopsomp' is often both the fastest and most accurate.<br>
<br>
James Tursa

Sat, 06 Nov 2010 00:25:05 +0000
Re: LARGE roundoff error in matrixvector multiplication?
http://www.mathworks.com/matlabcentral/newsreader/view_thread/295702#793633
Roger Stafford
"Henry Wolkowicz" <hwolkowicz@uwaterloo.ca> wrote in message <ib1q18$51r$1@fred.mathworks.com>...<br>
> I get the following relatively large roundoff error on matrixvector multiplication<br>
> <br>
> >> [size(A) nnz(A) nnz(x) normest(round(A*1e3)/1e3A)]<br>
> ans =<br>
> 129 185 465 185 0<br>
> >> norm(A(:,colPerm)*x(colPerm)  A*x)<br>
> ans =<br>
> 1.3523e012<br>
> <br>
> i.e. a A is sparse 129by185 and I am finding A*x<br>
> 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?<br>
> Am I doing something wrong here?<br>
         <br>
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 53bit 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.<br>
<br>
Try this experiment.<br>
<br>
x = 54321*pi;<br>
y = x*(1/3)  x*(1/5)  x*(1/7) + x*(1/105)<br>
<br>
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.<br>
<br>
Roger Stafford

Sat, 06 Nov 2010 07:42:03 +0000
Re: LARGE roundoff error in matrixvector multiplication?
http://www.mathworks.com/matlabcentral/newsreader/view_thread/295702#793645
Bruno Luong
Please see an old thread about sum, order, and multithreading:<br>
<a href="http://www.mathworks.com/matlabcentral/newsreader/view_thread/260828">http://www.mathworks.com/matlabcentral/newsreader/view_thread/260828</a><br>
<br>
In here a serious bug with V 2009A bug has been reported. The bug has been fixed since 2009B.<br>
<a href="http://www.mathworks.com/matlabcentral/newsreader/view_thread/261541">http://www.mathworks.com/matlabcentral/newsreader/view_thread/261541</a><br>
<br>
Bruno