0

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.

Community
  • 1
  • 1
Ramin
  • 2,133
  • 14
  • 22

2 Answers2

2

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.

Acorbe
  • 8,367
  • 5
  • 37
  • 66
  • 1
    +1, but not the answer, to use Choleskey decomposition, **B** must be positive definite. – Ramin Dec 05 '12 at 18:31
  • 1
    but it's a calculation of eigenvalues, you have big matrix thus you loose accuracy. Enforce the symmetry by doing C = .5 * (C + C'), Is B 1e-16 symmetric? – Acorbe Dec 05 '12 at 18:47
1

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,--

laaph
  • 301
  • 1
  • 4
  • Thanks; I get 7.105427357601e-015 too; The elements of my matrices are larger than **a** and **b** in this example. – Ramin Dec 06 '12 at 14:59