possible bug in solve_right() with arbitrary precision real numbers

[Daniel Friedan]

solve_right() fails on matrices of arbitrary precision real numbers for some but not all values of the precision.  See examples below.

I'm seeing this in Sage 5.5 64bit under OS X 10.6.8 and also in Sage 5.2 under ScientificLinux.

sage: RF= RealField(52)
sage: Y=matrix(RF,2,2,1)
sage: A = matrix(RF,[[0.24,1,0],[1,0,0]])
sage: A.solve_right(Y)
    ValueError: matrix equation has no solutions

sage: RF= RealField(53)
sage: Y=matrix(RF,2,2,1)
sage: A = matrix(RF,[[0.24,1,0],[1,0,0]])
sage: A.solve_right(Y)

[ 0.000000000000000   1.00000000000000]
[  1.00000000000000 -0.240000000000000]
[ 0.000000000000000  0.000000000000000]

 

solve_right() fails for other values of the precision, including 10, 58, and 98.  I haven't tried to search systematically.

[Dima Pasechnik]

On 2013-01-09, Daniel Friedan <dfri...@gmail.com> wrote:
> ------=_Part_1296_8764740.1357697544661
> Content-Type: text/plain; charset=ISO-8859-1

>
> solve_right() fails on matrices of arbitrary precision real numbers for
> some but not all values of the precision. See examples below.
>
> I'm seeing this in Sage 5.5 64bit under OS X 10.6.8 and also in Sage 5.2
> under ScientificLinux.
>
> sage: RF= RealField(52)
> sage: Y=matrix(RF,2,2,1)
> sage: A = matrix(RF,[[0.24,1,0],[1,0,0]])
> sage: A.solve_right(Y)
> ValueError: matrix equation has no solutions

indeed, also fails for me (Sage 5.6.beta2 on OSX 10.6.8).
It's certainly a bug.

[Volker Braun]

By default, solve_right() checks that the result is true. This fails due to numerical errors depending on the precision:

sage: RF= RealField(52)

sage: Y=matrix(RF,2,2,1)

sage: A = matrix(RF,[[0.24,1,0],[1,0,0]])

sage: sol = A.solve_right(Y, check=False)

sage: A*sol

[     1.00000000000000 -5.55111512312578e-17]

[    0.000000000000000      1.000000000000000]

[Volker Braun]

This is now http://trac.sagemath.org/13932

[Daniel Friedan]

Should solve_right() be introducing a numerical error here?  The only numbers involved are 0, 1, and 0.24 and the calculation is linear in 0.24.

sage: RF= RealField(52)
sage: Y=matrix(RF,2,2,1);Y

[ 1.00000000000000 0.000000000000000]
[0.000000000000000  1.00000000000000]

sage: A = matrix(RF,[[0.24,1,0],[1,0,0]]); A

[0.240000000000000  1.00000000000000 0.000000000000000]
[ 1.00000000000000 0.000000000000000 0.000000000000000]

sage: X=A.solve_right(Y, check=false); X

[ 0.000000000000000   1.00000000000000]
[  1.00000000000000 -0.240000000000000]
[ 0.000000000000000  0.000000000000000]

sage: A*X-Y

[    0.000000000000000 -5.55111512312578e-17]
[    0.000000000000000     0.000000000000000]

[Volker Braun]

For starters, 0.24 can't be represented exactly. An oldie but goodie on the subject:

http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html

[Daniel Friedan]

I think I understand now.  The calculation is surely not linear in 0.24.  What appears in the solution is 1/(1/0.24) which is not exactly 0.24.

thanks

[Dima Pasechnik]

On 2013-01-09, Daniel Friedan <dfri...@gmail.com> wrote:

> ------=_Part_186_18436191.1357759234951
> Content-Type: text/plain; charset=ISO-8859-1

>
> I think I understand now. The calculation is surely not linear in 0.24.
> What appears in the solution is 1/(1/0.24) which is not exactly 0.24.

yes, those are the joys of floating point numbers
(and someone blindly using '==' on floating point numbers)

IMHO one might like a setting when == on such data throws an error.
At least for debugging.
This is actually doable with Sage...

Dmitrii