Thread Subject:

Rounding error amplification in matrix/matrix multiplication

Hello,
I am currently working on my own implementation in C++ of Strassen's fast matrix multiplication algorithm. It works actually pretty well.
What I do now is unit testing (cppunit), where I use Matlab as a golden standard. Everything works fine, except for big matrices full of big float numbers. Neither my new algorithm nor the former (O(n^3)) can stand the comparison with the matlab's result, the delta is just way too big for every and each element. I made some research, and it occured to me that I probably pull a rounding error all over my computation (rounding error amplification), and I made the assumption that somehow, Matlab's algorithm compensates this rounding error amplification.
I would be very happy if somebody here could provide some information about that.

Thank you,

Stéphane

"Stéphane " <jalisastef@yahoo.ca> wrote in message <hpnmqc$gu7$1@fred.mathworks.com>...
> Hello,
> Matlab's algorithm compensates this rounding error amplification.

Not that I'm aware of (in the strict sense). However Matlab uses BLAS for full matrix product, and there is one thing you should pay attention is the so called "floating point control word" (FPCW) that control the precision on the internal computation inside the (Intel) processor. There is an compilation option on MSVC to allow to switch the FPCW between two states, otherwise you can control finer within the program. I believe Matlab by default uses 53-bit internal calculation. You should at least use as much.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hpnorh$idd$1@fred.mathworks.com>...
> "Stéphane " <jalisastef@yahoo.ca> wrote in message <hpnmqc$gu7$1@fred.mathworks.com>...
> > Hello,
> > Matlab's algorithm compensates this rounding error amplification.
>
> Not that I'm aware of (in the strict sense). However Matlab uses BLAS for full matrix product, and there is one thing you should pay attention is the so called "floating point control word" (FPCW) that control the precision on the internal computation inside the (Intel) processor. There is an compilation option on MSVC to allow to switch the FPCW between two states, otherwise you can control finer within the program. I believe Matlab by default uses 53-bit internal calculation. You should at least use as much.
>
> Bruno

One additional note for OP. The BLAS that MATLAB uses was not written by The Mathworks but is actually a 3rd party library. You *might* be able to find out some info about it on the web somewhere but I am guessing that most of the details are proprietary.

James Tursa

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hpnorh$idd$1@fred.mathworks.com>...
> "Stéphane " <jalisastef@yahoo.ca> wrote in message <hpnmqc$gu7$1@fred.mathworks.com>...
> > Hello,
> > Matlab's algorithm compensates this rounding error amplification.
>
> Not that I'm aware of (in the strict sense). However Matlab uses BLAS for full matrix product, and there is one thing you should pay attention is the so called "floating point control word" (FPCW) that control the precision on the internal computation inside the (Intel) processor. There is an compilation option on MSVC to allow to switch the FPCW between two states, otherwise you can control finer within the program. I believe Matlab by default uses 53-bit internal calculation. You should at least use as much.
>
> Bruno

Hello Bruno,

you're right about the 53 bits... it is the same for MSVC however, the double format has a 53 bits mantissa.
Thank for the clue about MSVC, I changed the precision from precise to strict in the floating point optimization option. I will let you know about it, tests are under process right now.

Stéphane

"Stéphane " <jalisastef@yahoo.ca> wrote in message <hpnmqc$gu7$1@fred.mathworks.com>...
>
> I am currently working on my own implementation in C++ of Strassen's fast matrix multiplication algorithm.
> Everything works fine, except for big matrices full of big float numbers. Neither my new algorithm nor the former (O(n^3)) can stand the comparison with the matlab's result, the delta is just way too big for every and each element.

Can you provide a short code sample that generates the matrices that are giving you problems?

James Tursa

"Stéphane " <jalisastef@yahoo.ca> wrote in message <hpnpph$6ii$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hpnorh$idd$1@fred.mathworks.com>...
> > "Stéphane " <jalisastef@yahoo.ca> wrote in message <hpnmqc$gu7$1@fred.mathworks.com>...
> > > Hello,
> > > Matlab's algorithm compensates this rounding error amplification.
> >
> > Not that I'm aware of (in the strict sense). However Matlab uses BLAS for full matrix product, and there is one thing you should pay attention is the so called "floating point control word" (FPCW) that control the precision on the internal computation inside the (Intel) processor. There is an compilation option on MSVC to allow to switch the FPCW between two states, otherwise you can control finer within the program. I believe Matlab by default uses 53-bit internal calculation. You should at least use as much.
> >
> > Bruno
>
> Hello Bruno,
>
> you're right about the 53 bits... it is the same for MSVC however, the double format has a 53 bits mantissa.
> Thank for the clue about MSVC, I changed the precision from precise to strict in the floating point optimization option. I will let you know about it, tests are under process right now.
>
> Stéphane

Hello James,

I am not asking for code, all I want to know is if matlab's algorithm compensates rounding errors amplifications or not. I could not find any clue about it on the web.
Thanks for the answer though. :)

Stéphane

"James Tursa" <aclassyguy_with_a_k_not_a_c@hotmail.com> wrote in message <hpnpta$8ca$1@fred.mathworks.com>...
> "Stéphane " <jalisastef@yahoo.ca> wrote in message <hpnmqc$gu7$1@fred.mathworks.com>...
> >
> > I am currently working on my own implementation in C++ of Strassen's fast matrix multiplication algorithm.
> > Everything works fine, except for big matrices full of big float numbers. Neither my new algorithm nor the former (O(n^3)) can stand the comparison with the matlab's result, the delta is just way too big for every and each element.
>
> Can you provide a short code sample that generates the matrices that are giving you problems?
>
> James Tursa


Here is the matlab code:

temp1 = rand(2000, 2000);
temp2 = rand(2000, 2000);

A = 25000 * temp1;
B = 25000 * temp2;

C = A * B;

fichier = fopen('A.txt', 'w+');
fprintf(fichier, '%16f;', A);
fclose(fichier);

fichier = fopen('B.txt', 'w+');
fprintf(fichier, '%16f;', B);
fclose(fichier);

fichier = fopen('C.txt', 'w+');
fprintf(fichier, '%16f;', C);
fclose(fichier);

Then, I launch my program on these matrices: A and B to perform the multiplication, which is then compared to result C.


Stéphane

"Stéphane " <jalisastef@yahoo.ca> wrote in message <hpnqad$ebj$1@fred.mathworks.com>...
>
> >
> > you're right about the 53 bits... it is the same for MSVC however, the double format has a 53 bits mantissa.

Sthéphane,

To avoid any confusion, the FPCW is an independent control on top of the size of the mantissa of the floating points. It controls how the processor perform internally the round off when operating on the floating point. For example, you could force the internal calculation to be 64-bits on 53-bit-mantissa numbers.

Bruno

Take a look at this link if you are not familiar with the FPCW:

http://msdn.microsoft.com/en-us/library/aa289157%28VS.71%29.aspx

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hpnra9$s8r$1@fred.mathworks.com>...
> Take a look at this link if you are not familiar with the FPCW:
>
> http://msdn.microsoft.com/en-us/library/aa289157%28VS.71%29.aspx
>
> Bruno

Yes I was on this page, that's why I told you about fp:fast, fp:strict and fp:precise (strict and precise :) ) in a previous message :)
No change though... :(

thank you,

Stéphane

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <hpnra9$s8r$1@fred.mathworks.com>...
> Take a look at this link if you are not familiar with the FPCW:
>
> http://msdn.microsoft.com/en-us/library/aa289157%28VS.71%29.aspx
>
> Bruno

Yes I was on this page, that's why I told you about fp:fast, fp:strict and fp:precise (strict and precise :) ) in a previous message :)
No change though... :(

thank you,

Stéphane

Dear Stéphane!

> I am not asking for code, all I want to know is if matlab's algorithm compensates rounding errors amplifications or not. I could not find any clue about it on the web.

Matlab uses BLAS for matrix (and array) operations and in consequence no error compensation is applied.
You can compile XBLAS and tell Matlab to use this instead. XBLAS uses the TwoSum, FastTwoSum and Split algorithms to calculate error compensated dot products and matrix multiplications. The results of the dot products are as accurate as if they are calculated with quadruple precision (but rounded to doubles finally of course).

The BLAS is compiled such that it uses FMAC commands: Fused Multiply Accumulate, which means that t = r + s * q is calculated with a single command and a single rounging error only, see:
  http://www.cs.berkeley.edu/~wkahan/MxMulEps.pdf
Kahan states:
"A less obvious advantage of a FMAC is slightly more accurate matrix products, which tend to an extra sig. bit of accuracy because they accumulate about half as many rounding errors."

A Matlab implementation of error compensated array functions is Rump's INTLAB: http://www.ti3.tu-harburg.de/rump/intlab/

Kind regards, Jan

"Jan Simon" <matlab.THIS_YEAR@nMINUSsimon.de> wrote in message <hpnsr1$mnt$1@fred.mathworks.com>...
> Dear Stéphane!
>
> > I am not asking for code, all I want to know is if matlab's algorithm compensates rounding errors amplifications or not. I could not find any clue about it on the web.
>
> Matlab uses BLAS for matrix (and array) operations and in consequence no error compensation is applied.
> You can compile XBLAS and tell Matlab to use this instead. XBLAS uses the TwoSum, FastTwoSum and Split algorithms to calculate error compensated dot products and matrix multiplications. The results of the dot products are as accurate as if they are calculated with quadruple precision (but rounded to doubles finally of course).
>
> The BLAS is compiled such that it uses FMAC commands: Fused Multiply Accumulate, which means that t = r + s * q is calculated with a single command and a single rounging error only, see:
> http://www.cs.berkeley.edu/~wkahan/MxMulEps.pdf
> Kahan states:
> "A less obvious advantage of a FMAC is slightly more accurate matrix products, which tend to an extra sig. bit of accuracy because they accumulate about half as many rounding errors."
>
> A Matlab implementation of error compensated array functions is Rump's INTLAB: http://www.ti3.tu-harburg.de/rump/intlab/
>
> Kind regards, Jan

Thank you very much Jan,

but, if there are no compensations applied, why are the results so differents? I mean, if it was only my implementation of Strassen, but the ancient O(n^3) 3-for-loops-matrix-multiplication-algorithm is far from Matlab too... (but very close to my Strassen...).
I do not know where I could find an answer now. The compensation was my last idea.
It is clearly a rounding error amplification problem though, since for small matrices I am able to achieve the desired precision, in the same test conditions.

Stéphane

Dear Stéphane!

This is in interesting paper of Kahan also:
  http://www.cs.berkeley.edu/~wkahan/Qdrtcs.pdf

He uses a small matrix operations to identify the rounding behaviour of Matlab:
  x = [1+4.656612873e-10, 1] * [1-4.656612873e-10; -1] ;
  y = [1, 1+4.656612873e-10] * [-1; 1-4.656612873e-10] ;

"If x and y are both zero, MATLAB follows the first style, rounding every arithmetic operation to the same 8-byte precision as is normally used for variables stored in memory. If x and y are both –2–62 » –2.1684e–19 , MATLAB follows the second style when computing scalar products of vectors and of rows and columns of matrices; these are accumulated to 64 sig. bits in 10-byte registers before being rounded back to 53 sig. bits when stored into 8-byte variables.

MATLAB 6.5 on PCs follows the first style by default; to enable the second style execute the command
  system_dependent( 'setprecision', 64 )
A similar command
  system_dependent( 'setprecision', 53 )
restores the default first style."

You can test the sensitivity of your matrix computations by comparing
  A * B and (A+eps) * B

Good luck, Jan

"Stéphane " <jalisastef@yahoo.ca> wrote in message <hpnqll$j8d$1@fred.mathworks.com>...
>
> Here is the matlab code:
>
> temp1 = rand(2000, 2000);
> temp2 = rand(2000, 2000);
>
> A = 25000 * temp1;
> B = 25000 * temp2;
>
> C = A * B;
>
> fichier = fopen('A.txt', 'w+');
> fprintf(fichier, '%16f;', A);
> fclose(fichier);
>
> fichier = fopen('B.txt', 'w+');
> fprintf(fichier, '%16f;', B);
> fclose(fichier);
>
> fichier = fopen('C.txt', 'w+');
> fprintf(fichier, '%16f;', C);
> fclose(fichier);
>
> Then, I launch my program on these matrices: A and B to perform the multiplication, which is then compared to result C.

For random matrices like these I wouldn't expect anything unusual in the accuracy of the result. I spot checked the MATLAB result against a high precision routine I have and saw the expected relative difference between the two results in the e-16 range. However, to be sure you are comparing apples to apples why don't you write and read a binary file instead of a text file? It may be pointless to talk about rounding modes of the processor if your data transfer via a text file write is introducing rounding differences in the data right up front, and I suspect this is where your real differences are coming from. What kind of relative differences are you seeing in your results vs MATLAB? Also, the 25000 factor isn't really doing anything for this test ... all of the relative sizes of the individual products are the same. Multiplying by 25000 simply has the qualitative effect of changing
all the exponents of all the values by the same amount, but since the relative sizes remain the same between them the accuracy of the matrix product will be the same.

James Tursa

"Stéphane " <jalisastef@yahoo.ca> wrote in message <hpnmqc$gu7$1@fred.mathworks.com>...
> Hello,
> I am currently working on my own implementation in C++ of Strassen's fast matrix multiplication algorithm. It works actually pretty well.
> What I do now is unit testing (cppunit), where I use Matlab as a golden standard. Everything works fine, except for big matrices full of big float numbers. Neither my new algorithm nor the former (O(n^3)) can stand the comparison with the matlab's result, the delta is just way too big for every and each element. I made some research, and it occured to me that I probably pull a rounding error all over my computation (rounding error amplification), and I made the assumption that somehow, Matlab's algorithm compensates this rounding error amplification.
> I would be very happy if somebody here could provide some information about that.
>
> Thank you,
>
> Stéphane

-------------------------------------

Dear Stéphane,

I recommend reading Chapter 22 - Fast Matrix Multiplication, Nicholas J Higham, "Accuracy and Stability of Numerical Algorithms", 2nd Edition, SIAM, 2002.

At the top of page 442 he gives this example :

A = [1, 0; 0, 1]; B = [1, e; e, e^2].

He says that C = A*B is calculated exactly in floating point by the conventional algorithm, whereas Strassen's algorithm computes

c22 = 2(1+e^2)+(e-e^2)-1-(1+e)

which has a relative error O(u/e^2), (u =eps/2) and is much larger than u if e is small.

He then gives this example, X = P32*E, where Pn is an n x n Pascal Matrix (in his Matrix Toolbox) and E = [eij] = [1/3]. He says that with just one level of recursion in Strassen's method, the max rel. err. of xij is about 10^(-5).

Derek O'Connor

"Derek O'Connor" <derekroconnor@eircom.net> wrote in message <hpoj5s$l0k$1@fred.mathworks.com>...
> "Stéphane " <jalisastef@yahoo.ca> wrote in message <hpnmqc$gu7$1@fred.mathworks.com>...
> > Hello,
> > I am currently working on my own implementation in C++ of Strassen's fast matrix multiplication algorithm. It works actually pretty well.
> > What I do now is unit testing (cppunit), where I use Matlab as a golden standard. Everything works fine, except for big matrices full of big float numbers. Neither my new algorithm nor the former (O(n^3)) can stand the comparison with the matlab's result, the delta is just way too big for every and each element. I made some research, and it occured to me that I probably pull a rounding error all over my computation (rounding error amplification), and I made the assumption that somehow, Matlab's algorithm compensates this rounding error amplification.
> > I would be very happy if somebody here could provide some information about that.
> >
> > Thank you,
> >
> > Stéphane
>
> -------------------------------------
>
> Dear Stéphane,
>
> I recommend reading Chapter 22 - Fast Matrix Multiplication, Nicholas J Higham, "Accuracy and Stability of Numerical Algorithms", 2nd Edition, SIAM, 2002.
>
> At the top of page 442 he gives this example :
>
> A = [1, 0; 0, 1]; B = [1, e; e, e^2].
>
> He says that C = A*B is calculated exactly in floating point by the conventional algorithm, whereas Strassen's algorithm computes
>
> c22 = 2(1+e^2)+(e-e^2)-1-(1+e)
>
> which has a relative error O(u/e^2), (u =eps/2) and is much larger than u if e is small.
>
> He then gives this example, X = P32*E, where Pn is an n x n Pascal Matrix (in his Matrix Toolbox) and E = [eij] = [1/3]. He says that with just one level of recursion in Strassen's method, the max rel. err. of xij is about 10^(-5).
>
> Derek O'Connor

Thank you everybody!

Your answers are very helpful, however I guess it will take me a few days to get back to you... :)

"James Tursa" <aclassyguy_with_a_k_not_a_c@hotmail.com> wrote in message <hpo19j$6n3$1@fred.mathworks.com>...
>
> It may be pointless to talk about rounding modes of the processor if your data transfer via a text file write is introducing rounding differences in the data right up front, and I suspect this is where your real differences are coming from.

e.g., try this using your format:

>> fid = fopen('A.txt', 'w+');
>> fprintf(fid, '%16f;', A); % Using your format
>> fclose(fid);
>> fid = fopen('A.txt', 'r');
>> B = fscanf(fid,'%16f;',inf);
>> max(abs(A(:)-B(:)))
ans =
  5.0000e-007

You are not starting with the same data for your comparison, so it is not a valid comparison.

Now try this, doing the writing and reading in binary:

>> fid = fopen('A.bin', 'w');
>> fwrite(fid,A,'double');
>> fclose(fid);
>> fid = fopen('A.bin', 'r');
>> B = fread(fid,inf,'double');
>> max(abs(A(:)-B(:)))
ans =
     0

Here you have the same data for comparison.

James Tursa

"Jan Simon" <matlab.THIS_YEAR@nMINUSsimon.de> wrote in message <hpo10m$27d$1@fred.mathworks.com>...
> He uses a small matrix operations to identify the rounding behaviour of Matlab:
> x = [1+4.656612873e-10, 1] * [1-4.656612873e-10; -1] ;
> y = [1, 1+4.656612873e-10] * [-1; 1-4.656612873e-10] ;
>
> "If x and y are both zero, MATLAB follows the first style, rounding every arithmetic operation to the same 8-byte precision as is normally used for variables stored in memory. If x and y are both –2–62 » –2.1684e–19 , MATLAB follows the second style when computing scalar products of vectors and of rows and columns of matrices; these are accumulated to 64 sig. bits in 10-byte registers before being rounded back to 53 sig. bits when stored into 8-byte variables.
>
> MATLAB 6.5 on PCs follows the first style by default; to enable the second style execute the command
> system_dependent( 'setprecision', 64 )
> A similar command
> system_dependent( 'setprecision', 53 )
> restores the default first style."
>
> You can test the sensitivity of your matrix computations by comparing
> A * B and (A+eps) * B
>
> Good luck, Jan

I'm calculating the norm of an array of size 3. I'm running my code on Matlab 2010a 32-bit and 64-bit on the same physical machine, where I'm getting x=0, y=0. On the 32-bit matlab, the result is one bit appart from the result on the 64-bit version. This small discrepancy adds up when the code is in a long loop. On a Windows environment 32-bit and 64-bit I get yet other results. The difference is likely coming from the order the actual operations are performed. Is there a way to make MATLAB behave the same on all these different environments?

Here is my sample code, using "format hex":
a = [hex2num('3fa6d7a56de3326b') hex2num('bfb23ec02f2f9874') hex2num('3f92630ec311622c')]
normA = norm(a)

32-bit:
a =
   3fa6d7a56de3326b bfb23ec02f2f9874 3f92630ec311622c
normA =
   3fb602a9c1574f01

64-bit:
a =
   3fa6d7a56de3326b bfb23ec02f2f9874 3f92630ec311622c
normA =
   3fb602a9c1574f02