# Is order of matrix multiplication affecting numerical accuracy of the result?

I have to multiply three matrices of floats: A (100x8000), B (8000x27) and C (27x1).

Is there any difference in accuracy between A(BC) and (AB)C? If yes - how may I determine the more accurate multiplication order? Speed is not a factor here.

Matrices A and B contain 8000 samples of (respectively) 100 and 27 features.

Yes, there are differences in accuracy since with machine numbers the usual properties of arithmetics don't hold.

Machine numbers are defined as $$F(\beta,t,m,M)= \{ 0 \} \cup \{ x \in \mathbb{R} : x = sign(x)\beta^{p} \sum_{i=1}^{t}d_i\beta^{-i},\ 0 \leq d_i \lt \beta\ ,\ d_1\ne 0\ , -m \le p \le M \}$$ and represent the subset of $$\mathbb{R}$$ that your machine is able to represent.

All other numbers must be approximated with a number in this subset (usually by truncating the numbers or rounding them).

Let's assume we are in $$F(10,2,m,M)$$ meaning we are working in base 10 with two digits.

Let $$x=0.11*10^1$$, $$y=0.31*10^{1}$$ and $$z=0.25*10^{1}$$.

The associative property of multiplication doesn't hold:

$$(x*y)*z = 0.34*10^1 * z = 0.85*10^1$$

and

$$x*(y*z) = x * (0.78 * 10^1) = 0.86*10^1$$

other properties that don't hold are:

• the associative property of the addition

• distributive properties

• $$x(y/x)$$ isn't always equals to $$x$$
• if $$xy=yz$$ then isn't always true that $$x=z$$
• Notice how $0.11 * 0.31 = 0.341$ and since $t = 2$ the $1$ is truncated. The operation on machine numbers are also called machine operations – Gerardo Zinno Jul 8 at 18:41
• For additional details: "What every computer scientist should know about floating-point arithmetic". The upshot is that getting floating-point arithmetic right is challenging, as Gerardo Zinno's answer illustrates. – user888379 Jul 8 at 19:36
• But is there any rule in which order the matrix multiplication should be performed for maximum accuracy? I can prove that for matrices which all elements have the same sign the error estimate is not affected by the order of multiplication. I have no idea whether that holds for any matrices (due to the possible loss of significant digits). – abukaj Jul 9 at 11:33
• @abukaj I'm not sure, but I think there isn't a way to tell a priori which way is gonna be more precise since it depends on which numbers are in the matrix. The only thing you could in theory do is multiply the matrices by hand or on a system with greater precision and then compare the result to see which way was better. But then again, with different matrices the you would have to repeat it all again. Also, you would die before you could multiply by hand those matrices. – Gerardo Zinno Jul 9 at 13:00
• You could also use a system with arbitrary precision (given sufficient memory), instead of a system with greater precision. I have Mathematica in mind. – Simon Jul 12 at 10:49

Interesting problem. Obviously with floating point involved, you can’t expect the same results, but we’d like to know what tended to give better results?

In a completely unscientific way, you can perform both methods with extended precision and double precision and compare the results. You’d hope that the extended precision results are much closer together, and would guess that the double precision result closer to both is the better one.

If the numbers are small, you could round each to a multiple of 0.1, calculate the product, and since each item in the result must be a multiple of 0.001 you might be able to figure out the exact rounding error.

At the moment the best solution I have (for a particular set of matrices) is to perform the following numerical experiment:

1. Calculate a reference matrix as an average of products calculated with high precision (e.g. `np.float128).
2. Calculate test products with lower precision (np.float64, np.float32, even np.float16),
3. Analyse errors calculated as a difference between test products and the reference matrix. The errors are expected to decline as the precision is higher.
• I may be missing something, but it seems that this is only a proof that you are able to get the same error bound on both methods. This may be the best you can do in your model, but there may be other effects not included in you model that lead to significant difference in practice. – Discrete lizard Jul 11 at 15:26
• @Discretelizard You are right, there were serious issues with the proof. – abukaj Jul 11 at 15:47

how may I determine the more accurate multiplication order?

This is a bit of a frame challenge. To determine the more accurate order you can perform both orders with interval arithmetic. But having done that, you can intersect the intervals to get a more precise answer than either individual order.

Of course, the dot product used will probably be a more important factor than the order of multiplication.