CODEWITHSUNDEEP
matlabmatrix-multiplication
Ramin
Ramin

Reputation: 2133

Accurate matrix multiplication in Matlab

There are 2 matrices:

A: (6 x 78) max=22.2953324329113, min=0
B: (6 x 6 ) max=2187.9013214004 , min=-377.886378385521

B is symmetric and as a result, C = A' * B * A must be a symmetric matrix (theoretically), but this is not the case when I calculate them in Matlab. In fact:

max(max(abs(C - C'))) = 2.3283064365386963e-010

How can I multiply them and get an accurate result?
or
What is a safe way to round the elements of C?

I read this question : efficient-multiplication-of-very-large-matrices-in-matlab, but my problem is not speed or memory. I need an accurate result

Thanks.

Upvotes: 0

Views: 1659

Answers (2)

laaph
laaph

Reputation: 301

You need to read "What Every Computer Scientist Should Know About Floating-Point Arithmetic":

http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html

Realize that computers will never be able to give perfect floating point results, and that leaves you with a few options:

  • Do as few operations as possible, that is, pick the order of operations for the purpose of having the fewest rounding errors
  • Fixed decimal point arithmetic - or integer arithmetic - this isn't always practical for all applications, but in some applications, you can get away with this. Financial applications are the commonly cited example (multiply by 100 to make pennies go away! divide by 100 when you are done!).
  • There are other tricks I can't think of this late.

I'm going to have to give your operations a spin - on my machine, eps gives me 2.2204e-16, which is six orders of magnitude lower than what you are getting. See what eps is on your machine - it should be similar - if it is something like 1e-12 or so, I'd say your result is exactly what you'd expect from those operations.

When I do this with random numbers, I get

a = rand(6, 78);
b = rand(6, 6);
b = b + b';   % To make b symmetric
c = a' * b * a;
max(max(abs(c - c')))

ans =

   7.1054e-15

Which is a bit closer to what I'd expect with rounding errors after that many operations, but I am not sure of your input, your machine, and I have no idea what else might be affecting things.

Cheers,--

Upvotes: 1

Acorbe
Acorbe

Reputation: 8391

You can consider cholesky decomposition of B since it is symmetric

  B = R'R
  R = chol(A)   % // in matlab

then C = A'R'R A =D'D, where D = RA.

With C=D'D you should have machine epsilon precision, although you introduce a possible error due to the accuracy of the decomposition.

Upvotes: 2

Related Questions