3
$\begingroup$

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)$$

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Federico Poloni Sep 22 at 12:33
3
$\begingroup$

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.

Workarounds: either higher-precision or integer arithmetic, or in theory since you are computing with integers you could perform the computation modulo a few different coprime values and use the Chinese remainder theorem to infer the final result --- but if your computations are at this scale probably it's just not worth the effort.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.