Prev: using mex/fftw for "xcorr", faster than 2010A "xcorr"
Next: Fitting ellipses with dispersion
From: Derek O'Connor on 9 Apr 2010 21:10
"Stéphane " <jalisastef(a)yahoo.ca> wrote in message <hpnmqc$gu7$1(a)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
From: Stéphane on 12 Apr 2010 10:10
"Derek O'Connor" <derekroconnor(a)eircom.net> wrote in message <hpoj5s$l0k$1(a)fred.mathworks.com>...
> "Stéphane " <jalisastef(a)yahoo.ca> wrote in message <hpnmqc$gu7$1(a)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... :)
From: James Tursa on 12 Apr 2010 11:23
"James Tursa" <aclassyguy_with_a_k_not_a_c(a)hotmail.com> wrote in message <hpo19j$6n3$1(a)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
First  |  Prev  | 
Pages: 1 2 3 4
Prev: using mex/fftw for "xcorr", faster than 2010A "xcorr"
Next: Fitting ellipses with dispersion