Why does matrix power give wrong answer?

I am trying to do a MCMC simulation and need to calculate a small frequency event over a long time. To do this I need to raise

A=[1-1.01E-11,1.01E-11,0;1-3.024E-11,0,3.024E-11;0,0,1];

to a large power (10^10). I am particularly interested in the (1,3) element of the result (the probability of transferring between states 1 and 3 after 10 billion iterations). After 10^9 the (1,3) element is 3.05423999691491e-13. However after 10^10 it is -3.05423999991124e-22. Which is impossible because I am multiplying matrices with strictly non-negative values. As a check I tried (A^(10^9))^10, and this gave a (1,3) entry of 3.05423999966373e-12. Doing the calculation (A^(10^10) using repeated squaring and multiplication) in Excel gave 3.05424E-12.

I am using MATLAB 2013b on windows 7 64-bit.

Why does A^(10^10) fail but (A^(10^9))^10 work?

[John D'Errico]

Time to learn about floating point arithmetic. When your numbers are only held accurately to roughly 16 digits, what do you expect?

Answer 1

[per isakson]

See

• What Every Programmer Should Know About Floating-Point Arithmetic
• What Every Computer Scientist Should [...]
• Why is 0.3 - 0.2 - 0.1 (or similar) not equal to zero?
• and HPF - a big decimal class by John D'Errico

[Joseph]

Try

A=[1-1e-12,1e-12,0;1-1e-12,0,1e-12;0,0,1];
B=A^(2^31-1)*A

You get

0.999999999999000 9.99999999999000e-13 2.14748364699785e-15

0.999999999998000 9.99999999998000e-13 1.00214748364600e-12

0 0 1

Now try

B=A^(2^31)

You get

0.999999999999000 9.99999999999000e-13 -9.99999999997652e-25

0.999999999998000 9.99999999998000e-13 9.99999999997652e-13

0 0 1

The answers should be identical. As far as the floating point arithmetic argument, the second calculation can be done with 31 squarings, yet some other algorithm is used and the value in the (1,3) element is horribly wrong. Floating point arithmetic errors alone do not explain the difference in values between A^(2^31-1)*A and A^(2^31).

[per isakson]

"The answers should be identical." If you think Matlab does not behave as promised, you might want to report it as a bug.

[per isakson]

See MATLAB: Repeated multiplication is faster than integer powers together with its comments

[per isakson]

See Floating point addition is not associative