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:
LARGE roundoff error in matrix-vector multiplication?

Subject: LARGE roundoff error in matrix-vector multiplication?

From: Henry Wolkowicz

Date: 5 Nov, 2010 20:39:04

Message: 1 of 7

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?

Subject: LARGE roundoff error in matrix-vector multiplication?

From: Bruno Luong

Date: 5 Nov, 2010 22:16:03

Message: 2 of 7

"Henry Wolkowicz" <hwolkowicz@uwaterloo.ca> wrote in message <ib1q18$51r$1@fred.mathworks.com>...
> 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?

http://en.wikipedia.org/wiki/Condition_number

Your matrix is badly scaled.

Bruno

Subject: LARGE roundoff error in matrix-vector multiplication?

From: Henry Wolkowicz

Date: 5 Nov, 2010 23:04:04

Message: 3 of 7

Bruno ... your reply does not make sense to me. I am doing matrix-vector multiplication. What does an ill-conditioned matrix (if it is or not) have anything to do with matrix multiplication?
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.

 1.0e+003 *

   0.001000000000000 0.000000000000000 0.001000000000000 0.000000000000000
  -0.000400000000000 0.158000000000000 -0.001700000000000 0.158000000000000
  -0.000400000000000 0.158000000000000 0.001500000000000 2.504373333333334
   0.001500000000000 2.504373333333334 -0.000400000000000 0.158000000000000
  -0.001700000000000 0.158000000000000 -0.000400000000000 0.158000000000000
  -0.001700000000000 0.158000000000000 -0.001700000000000 0.158000000000000

when I multiply (inner product) the first two vectors and the last two I get

valrow3 =

    3.092960000000001e+003


valrow3perm =

    3.092960000000001e+003


ans =

    4.547473508864641e-013


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
of 1e-12




"Henry Wolkowicz" <hwolkowicz@uwaterloo.ca> wrote in message <ib1q18$51r$1@fred.mathworks.com>...
> 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?

Subject: LARGE roundoff error in matrix-vector multiplication?

From: Bruno Luong

Date: 5 Nov, 2010 23:26:04

Message: 4 of 7

"Henry Wolkowicz" <hwolkowicz@uwaterloo.ca> wrote in message <ib22h4$9bp$1@fred.mathworks.com>...
> Bruno ... your reply does not make sense to me. I am doing matrix-vector multiplication. What does an ill-conditioned matrix (if it is or not) have anything to do with matrix multiplication?

You are right, sorry somehow I though you solve linear system.


>
> valrow3 =
>
> 3.092960000000001e+003
>
>
> valrow3perm =
>
> 3.092960000000001e+003
>
>
> ans =
>
> 4.547473508864641e-013
>
>
> 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
> of 1e-12
>

4.547473508864641e-013 is eps(3.092960000000001e+003), so what you get is simply at a limit of floating-point accuracy. Is that surprise you? Not me.

Bruno

Subject: LARGE roundoff error in matrix-vector multiplication?

From: James Tursa

Date: 5 Nov, 2010 23:48:05

Message: 5 of 7

"Henry Wolkowicz" <hwolkowicz@uwaterloo.ca> wrote in message <ib22h4$9bp$1@fred.mathworks.com>...
> Bruno ... your reply does not make sense to me. I am doing matrix-vector multiplication. What does an ill-conditioned matrix (if it is or not) have anything to do with matrix multiplication?
> 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.
>
> 1.0e+003 *
>
> 0.001000000000000 0.000000000000000 0.001000000000000 0.000000000000000
> -0.000400000000000 0.158000000000000 -0.001700000000000 0.158000000000000
> -0.000400000000000 0.158000000000000 0.001500000000000 2.504373333333334
> 0.001500000000000 2.504373333333334 -0.000400000000000 0.158000000000000
> -0.001700000000000 0.158000000000000 -0.000400000000000 0.158000000000000
> -0.001700000000000 0.158000000000000 -0.001700000000000 0.158000000000000
>
> when I multiply (inner product) the first two vectors and the last two I get
>
> valrow3 =
>
> 3.092960000000001e+003
>
>
> valrow3perm =
>
> 3.092960000000001e+003
>
>
> ans =
>
> 4.547473508864641e-013
>
>
> 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
> of 1e-12

You can try these different methods to see numerical differences just by reordering the calculations:

Pick of two of your vectors for the dot product calculation, call them A and B, and then try these methods:

dot(A,B)
A.'*B
mtimesx(A,'t',B,'blas')
mtimesx(A,'t',B,'loops')
mtimesx(A,'t',B,'loopsomp')

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:

http://www.mathworks.com/matlabcentral/fileexchange/25977-mtimesx-fast-matrix-multiply-with-multi-dimensional-support

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.

James Tursa

Subject: LARGE roundoff error in matrix-vector multiplication?

From: Roger Stafford

Date: 6 Nov, 2010 00:25:05

Message: 6 of 7

"Henry Wolkowicz" <hwolkowicz@uwaterloo.ca> wrote in message <ib1q18$51r$1@fred.mathworks.com>...
> 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

Subject: LARGE roundoff error in matrix-vector multiplication?

From: Bruno Luong

Date: 6 Nov, 2010 07:42:03

Message: 7 of 7

Please see an old thread about sum, order, and multithreading:
http://www.mathworks.com/matlabcentral/newsreader/view_thread/260828

In here a serious bug with V 2009A bug has been reported. The bug has been fixed since 2009B.
http://www.mathworks.com/matlabcentral/newsreader/view_thread/261541

Bruno

Tags for this Thread

No tags are associated with 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