# Bounding error of float32 matrix multiplication

Some numerical debugging led me to the minimal example below. I'm observing relative error of 0.75 on individual elements.

Is there a way to estimate/bound this error without resorting to higher precision arithmetic?

import numpy as np
x = np.array([[11041, 13359, 15023, 18177], [13359, 16165, 18177, 21995], [15023, 18177, 20453, 24747], [18177, 21995, 24747, 29945]])
y = np.array([[29945, -24747, -21995, 18177], [-24747, 20453, 18177, -15023], [-21995, 18177, 16165, -13359], [18177, -15023, -13359, 11041]])

print(x@y)          # high precision

#[[16  0  0  0]
# [ 0 16  0  0]
# [ 0  0 16  0]
# [ 0  0  0 16]]

x0 = x.astype(np.float32)
y0 = y.astype(np.float32)
print(x0@y0)        # low precision

#[[28.  0. -4.  0.]
# [ 0. 28.  0. -4.]
# [20.  0. 20.  0.]
# [ 0. 20.  0. 20.]]



This example comes from using exact arithmetic to compute $$X^{-1}$$, and checking the value of $$XX^{-1}$$ in float32 where $$X = AA' \otimes BB'$$ with $$A=\left(\begin{matrix}5 & 6 \\ 7 & 8\end{matrix}\right)\ \ B=\left(\begin{matrix}9 & 10 \\ 11 & 12\end{matrix}\right)$$

• It's not clear to me what that notation with double letters means ($XX$, $AA$, $BB$), but it seems confusing at least. Also $\otimes$ is usually reserved for Kronecker product. Commented Sep 22, 2019 at 12:33

Estimating the error: yes, this is standard material in rounding error analysis. See for instance (3.13) in Higham's Accuracy and Stability of Numerical Algorithms: if $$C=AB$$ and $$\hat{C}$$ is its computed version, then $$|C-\hat{C}| \leq \gamma_n |A|\,|B|,$$ where the absolute values are taken componentwise and $$\gamma_n = \frac{nu}{1-nu} = nu + O(u^2)$$, with $$u$$ the machine precision (about 6e-8 in single precision) and $$n$$ the matrix size.
Essentially, the rounding error depends on the product of the absolute value matrices, or on the magnitude of intermediate results encountered when computing the scalar products, since you have to store those in single-precision variables. In your case, these intermediate values are of the order of about $$10^5 \cdot 10^5 = 10^{10}$$, so you can't expect smaller error than what you are getting. Note that a perturbation of your input values in their last (single-precision) significant digit would have perturbed the result by the same magnitude, so you can't really do better. The moment you write x0 = x.astype(np.float32) and forget that these are exact integers you have already lost.