![](https://secure.gravatar.com/avatar/744f2d519373c4944145986e434c0a79.jpg?s=120&d=mm&r=g)
Hi all, We are developing a code that heavily relies on NumPy. Some of our regression tests rely on floating point number comparisons and we are a bit lost in determining how to choose atol and rtol (we are trying to do all operations in double precision). We would like to set atol and rtol as low as possible but still have the tests pass if we run on different architectures or introduce some ‘cosmetic’ changes like using different similar NumPy routines. For example, let’s say we want some powers of the matrix A and compute them as: A = np.array(some_array) A2 = np.dot(A, A) A3 = np.dot(A2, A) A4 = np.dot(A3, A) If we alternatively computed A4 like: A4 = np.linalg.matrix_power(A, 4), we get different values in our final outputs because obviously the operations are not equivalent up to machine accuracy. Is there any reference that one could share providing guidelines on how to choose reasonable values for atol and rtol in this kind of situation? For example, does the NumPy package use a fixed set of values for its own development? the default ones? Thanks in advance for any help, Cheers, Bernard.
![](https://secure.gravatar.com/avatar/5f88830d19f9c83e2ddfd913496c5025.jpg?s=120&d=mm&r=g)
On Wed, Feb 24, 2021 at 11:29 AM Bernard Knaepen <bknaepen@gmail.com> wrote:
I don't think there's a clear guide in docs or blog post anywhere. You can get a sense of what works by browsing the unit tests for numpy and scipy. numpy.linalg, scipy.linalg and scipy.special are particularly relevant probably. For a rough rule of thumb: if you test on x86_64 and precision is on the order of 1e-13 to 1e-16, then set a relative tolerance 10 to 100 times higher to account for other hardware, BLAS implementations, etc. Cheers, Ralf
![](https://secure.gravatar.com/avatar/81e62cb212edf2a8402c842b120d9f31.jpg?s=120&d=mm&r=g)
Matrix powers are annoyingly tricky to keep under control due to the fact that things to explode or implode rather quickly. In fact the famous quote from Moler, Van Loan "Unfortunately, the roundoff errors in the mth power of a matrix, say B^m ,are usually small relative to ||B||^m rather than ||B^m||" So thanks for nothing smart people. Then there is a typical bound on rounding errors of matrix multiplication which says if C is the exact result and Ce the result of our operation then for some number k this expression holds |C - Ce| <= k * |A| * |B| From that, hoping that matrix power won't result with an error too far from the manual multiplication, it is a matter of selecting a sensible k for atol and rtol=0. I would go about this as (some arbitrary constant I am randomly throwing in)*(matrix size n)* np.finfo(dtype).eps*norm(A, 1)**k As an example, get a matrix and artificially bloat the (0,0) entry n = 100 A = np.random.rand(n, n) A += np.diag([10.]+[0.]*99) A4 = np.linalg.matrix_power(A, 4) AA = A @ A @ A @ A print('Max entry error', np.max(np.abs(AA-A4))) print('My atol value', 100*n*np.finfo(A.dtype).eps*np.linalg.norm(A, 1)*4) This accidentally makes it a tight bound but depending on how wildly your A varies or how the spectrum of A is structured you might need to change these constants. On Wed, Feb 24, 2021 at 1:53 PM Kevin Sheppard <kevin.k.sheppard@gmail.com> wrote:
![](https://secure.gravatar.com/avatar/5f88830d19f9c83e2ddfd913496c5025.jpg?s=120&d=mm&r=g)
On Wed, Feb 24, 2021 at 11:29 AM Bernard Knaepen <bknaepen@gmail.com> wrote:
I don't think there's a clear guide in docs or blog post anywhere. You can get a sense of what works by browsing the unit tests for numpy and scipy. numpy.linalg, scipy.linalg and scipy.special are particularly relevant probably. For a rough rule of thumb: if you test on x86_64 and precision is on the order of 1e-13 to 1e-16, then set a relative tolerance 10 to 100 times higher to account for other hardware, BLAS implementations, etc. Cheers, Ralf
![](https://secure.gravatar.com/avatar/81e62cb212edf2a8402c842b120d9f31.jpg?s=120&d=mm&r=g)
Matrix powers are annoyingly tricky to keep under control due to the fact that things to explode or implode rather quickly. In fact the famous quote from Moler, Van Loan "Unfortunately, the roundoff errors in the mth power of a matrix, say B^m ,are usually small relative to ||B||^m rather than ||B^m||" So thanks for nothing smart people. Then there is a typical bound on rounding errors of matrix multiplication which says if C is the exact result and Ce the result of our operation then for some number k this expression holds |C - Ce| <= k * |A| * |B| From that, hoping that matrix power won't result with an error too far from the manual multiplication, it is a matter of selecting a sensible k for atol and rtol=0. I would go about this as (some arbitrary constant I am randomly throwing in)*(matrix size n)* np.finfo(dtype).eps*norm(A, 1)**k As an example, get a matrix and artificially bloat the (0,0) entry n = 100 A = np.random.rand(n, n) A += np.diag([10.]+[0.]*99) A4 = np.linalg.matrix_power(A, 4) AA = A @ A @ A @ A print('Max entry error', np.max(np.abs(AA-A4))) print('My atol value', 100*n*np.finfo(A.dtype).eps*np.linalg.norm(A, 1)*4) This accidentally makes it a tight bound but depending on how wildly your A varies or how the spectrum of A is structured you might need to change these constants. On Wed, Feb 24, 2021 at 1:53 PM Kevin Sheppard <kevin.k.sheppard@gmail.com> wrote:
participants (4)
-
Bernard Knaepen
-
Ilhan Polat
-
Kevin Sheppard
-
Ralf Gommers