Is there a numerically optimal order of matrix multiplication?

Question

TL;DR: The question is about multiplication ACCURACY

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

Since matrices B and C are constant and A is variable, I prefer to calculate it as: ABC = np.dot(A, np.dot(B, C)). However I wonder, that it may be numerically worse (in terms of accuracy) than np.dot(np.dot(a, B), C).

What may be important: matrices A and B contain 8000 samples of (respectively) 100 and 27 correlated features.

Is there a numerically optimal (in terms of accuracy) order of the multiplication? If yes - how may I determine it?

Special Case

It may be assumed that both A and B matrices are nonnegative. Moreover:

C = np.linalg.solve(cov(B, k), X)

where X is a 27x1 matrix of 27 (possibly correlated) random variables of unknown distribution, cov = lambda X, k: np.dot(X.T, X) + k * np.eye(X.shape[1]), and k is a nonnegative constant minimizing the expression:

sum((X[i, 0] - np.dot(np.dot(B[:, [i]].T, drop(B, i)),
                      np.linalg.solve(cov(drop(B, i), k),
                                      np.delete(X, i, axis=0))) **2
    for i in range(27))

The drop() function is defined as lambda X, i: np.delete(X, i, axis=1).

Even More Special Case

It may be assumed that np.cov(B.T, B) is a covariance matrix of X, which follows multivariate Gaussian distribution.

Answer 1

At the moment the best idea 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.