Relative tolerance in French Sage book

[rbeezer]

A doctest from the French Sage book just barely fails in SageMathCloud (Intel Haswell, iirc).

CC: @zimmermann6 @ClementPernet

Component: documentation

Keywords: french book doctest tolerance

Author: Rob Beezer

Branch/Commit: bf5b521

Reviewer: Marc Mezzarobba

Issue created by migration from https://trac.sagemath.org/ticket/18593

[rbeezer]

Author: Rob Beezer

[rbeezer]

Branch: u/rbeezer/tolerance-french-book

[rbeezer]

Commit: bf5b521

[rbeezer]

comment:1

Increased the tolerance by an order of magnitude, so test now passes.

Paul - feel free to cc the relevant section author! Mostly experimenting making a contribution from with SageMathCloud, but thought it best to kill this one (repeated) failure. Thanks, Rob

---

New commits:

bf5b521

18593: Relative tolerance on linear system in French Sage book

[rbeezer]

comment:2

Fixed typo in ticket title. ;-)

[mezzarobba]

Reviewer: Marc Mezzarobba

[zimmermann6]

comment:4

just curious, I'd like to investigate more. Here is what I get with Sage 6.8.beta1 on my workstation:

sage: A = matrix(RDF, [[1,3,2],[1,4,2],[0,5,2],[1,3,2]]); A
[1.0 3.0 2.0]
[1.0 4.0 2.0]
[0.0 5.0 2.0]
[1.0 3.0 2.0]
sage: b = vector(RDF, [1,2,3,4]); b
(1.0, 2.0, 3.0, 4.0)
sage: transpose(A)
[1.0 1.0 0.0 1.0]
[3.0 4.0 5.0 3.0]
[2.0 2.0 2.0 2.0]
sage: R = transpose(A)*b; R
(7.0, 38.0, 20.0)
sage: Z = transpose(A)*A; Z
[ 3.0 10.0  6.0]
[10.0 59.0 30.0]
[ 6.0 30.0 16.0]
sage: Z.solve_right(R)
(-1.5000000000000049, -0.5000000000000011, 2.750000000000004)
sage: from sage.misc.citation import get_systems
sage: get_systems('Z.solve_right(R)')
['numpy', 'scipy']

What do you get on the machine where the test did fail?

Paul

[ClementPernet]

comment:5

One should not expect the numerical garbage to be reproducible when changing architectures: even if basic operations are compliant with the IEEE norm, the way they are scheduled (by the compiler or programmer) impacts the numerical noise.
Looking at the condition number (1222), we can assume that up to 3 digits can be lost. So the bug was indeed that the tolerance was probably set to a too low value.

I ran directly the LAPACK degesv implementation of GotoBLAS-2.11 on various architectures and got the following values:

x = [-1.5000000000000067 ,-0.50000000000000389 ,2.7500000000000098] for Intel Haswell, SandyBridge and Westmere, AMD K8,
x = [-1.4999999999999918 ,-0.49999999999999678 ,2.7499999999999911] for AMD Bulldozer
x = [-1.5000000000000135 ,-0.50000000000000855 ,2.7500000000000213] for AMD K10

[zimmermann6]

comment:6

indeed the -0.50000000000000855 middle value on AMD K10 gives a relative error of 1.53e-14, which is larger than the previous tolerance of 1e-14.

[rbeezer]

comment:7

Paul - Sorry to have been incomplete on this one. Indeed, the failure looks just like the AMD K10 "middle one."

Marc - thanks for the review!

File "src/sage/tests/french_book/linsolve_doctest.py", line 86, in sage.tests.french_book.linsolve_doctest
Failed example:
    Z.solve_right(R)  # rel tol 1e-14
Expected:
    (-1.5000000000000044, -0.5000000000000009, 2.750000000000003)
Got:
    (-1.5000000000000135, -0.5000000000000085, 2.7500000000000213)
Tolerance exceeded in 1 of 3:
    -0.5000000000000009 vs -0.5000000000000085, tolerance 2e-14 > 1e-14

[vbraun]

Changed branch from u/rbeezer/tolerance-french-book to bf5b521