Numerical noise with optimize.py on SPARC

[Dr. David Kirkby]

I get a doctest failure here on Solaris 10 SPARC with
sage/numerical/optimize.py. I can't find this is trac, but is anyone aware of it
being listed as a bug?


Does anyone know an accurate answer to this? I tried to do it in Mathematica,
using the FindMaximum[] command, but I can't see to get that to return anything
other than machine precision, despite adding options I would have thought would
have increased the accuracy

In[7]:= N[FindMaximum[8 E^(-x) Sin[x] -1,{x,0,8},AccuracyGoal->200,
PrecisionGoal->200],100]

Out[7]= {1.57918, {x -> 0.785398}}

In[8]:= Precision[%]

Out[8]= MachinePrecision

Anyway, here is the failure, on Solaris 10 SPARC, with gcc 4.5.0 (not the gcc in
Sage).


sage -t --long -force_lib devel/sage/sage/numerical/backends/coin_backend.pyx
[11.0 s]
sage -t --long -force_lib devel/sage/sage/numerical/optimize.py
**********************************************************************
File
"/export/home/drkirkby/sage-5.3.rc0/devel/sage-main/sage/numerical/optimize.py",
line 135:
sage: find_local_maximum(fast_float(8*e^(-x)*sin(x) - 1, x), 0, 8)
Expected:
(1.579175535558..., 0.78539817...)
Got:
(1.5791755355586758, 0.78539816111340355)
GLPK Simplex Optimizer, v4.44
6 rows, 3 columns, 8 non-zeros
Preprocessing...
2 rows, 2 columns, 4 non-zeros
Scaling...
A: min|aij| = 2.400e+01 max|aij| = 5.000e+01 ratio = 2.083e+00
GM: min|aij| = 8.128e-01 max|aij| = 1.230e+00 ratio = 1.514e+00
EQ: min|aij| = 6.606e-01 max|aij| = 1.000e+00 ratio = 1.514e+00
Constructing initial basis...
Size of triangular part = 2
* 0: obj = -5.100000000e+01 infeas = 0.000e+00 (0)
* 1: obj = -5.225000000e+01 infeas = 0.000e+00 (0)
OPTIMAL SOLUTION FOUND
**********************************************************************
1 items had failures:
1 of 7 in __main__.example_2
***Test Failed*** 1 failures.
For whitespace errors, see the file
/export/home/drkirkby/.sage/tmp/buzzard-2913/optimize_17468.py
[27.6 s]

--
Dr. David Kirkby Ph.D C.Eng MIET

[Dr. David Kirkby]

On 09/ 3/12 11:20 PM, Dr. David Kirkby wrote:
> I get a doctest failure here on Solaris 10 SPARC with
> sage/numerical/optimize.py. I can't find this is trac, but is anyone
> aware of it being listed as a bug?

No answers, so it is now

http://trac.sagemath.org/sage_trac/ticket/13427

> Does anyone know an accurate answer to this? I tried to do it in
> Mathematica, using the FindMaximum[] command, but I can't see to get
> that to return anything other than machine precision, despite adding
> options I would have thought would have increased the accuracy

Again no answers, but I managed to work out how to do it in Mathematica.

In[3]:= FindMaximum[8 E^(-x) Sin[x] -1,{x,0,8},WorkingPrecision->50]

Out[3]= {1.5791755355586755940189354764810356787023224274029,

> {x -> 0.78539816339744830961566077712689095303902649201882}}

In[4]:=

(This result was computed by Mathematica 7 running on OpenSolaris on a quad core
Xeon x86 CPU, though the working precision would have meant Mathematica did not
use the FPU)

OK, so now I have what probably is a fairly accurate answer.

> Anyway, here is the failure, on Solaris 10 SPARC, with gcc 4.5.0 (not
> the gcc in Sage).
>
>
> sage -t --long -force_lib
> devel/sage/sage/numerical/backends/coin_backend.pyx
> [11.0 s]
> sage -t --long -force_lib devel/sage/sage/numerical/optimize.py
> **********************************************************************
> File
> "/export/home/drkirkby/sage-5.3.rc0/devel/sage-main/sage/numerical/optimize.py",
> line 135:
> sage: find_local_maximum(fast_float(8*e^(-x)*sin(x) - 1, x), 0, 8)
> Expected:
> (1.579175535558..., 0.78539817...)
> Got:
> (1.5791755355586758, 0.78539816111340355)

Note the expected value of the Sage test "(1.579175535558..., 0.78539817...)" is
actually incorrect in the second instance as it should be 0.78539816 and not
0.78539817. The result returned on the SPARC processor is more accurate in this
particular instance, where the absolute error is -2.84x10^-9. So the test needs
changing to be

Expected:
(1.579175535558..., 0.7853981...)

to take account of the various errors when this test is run on different CPUs.

This is

Dave

[Dr. David Kirkby]

On 09/ 4/12 08:20 AM, Dr. David Kirkby wrote:
> On 09/ 3/12 11:20 PM, Dr. David Kirkby wrote:
>> I get a doctest failure here on Solaris 10 SPARC with
>> sage/numerical/optimize.py. I can't find this is trac, but is anyone
>> aware of it being listed as a bug?
>
> No answers, so it is now
>
> http://trac.sagemath.org/sage_trac/ticket/13427

I've attached a patch to

http://trac.sagemath.org/sage_trac/ticket/13427

which needs review. It solves the numerical noise issue, which for once is due
to the fact the SPARC is giving a more accurate result than the Intel chips!

The patch also reduces the range of the search, which was originally from x=0 to
x=8. It was pointed out on the trac ticket by Nils Bruin that there is another
maximum, at 9 Pi/4 (7.06858).

Dave