Disallow automatic abs() simplifications resulting non-real expressions

[jdemeyer]

Calling abs() on certain expressions can yield something which is not guaranteed to be evaluated as a real number:

sage: abs(x^2)
x*conjugate(x)
sage: abs(x)^2
x*conjugate(x)

In the presence of numerical noise, the expression x*conjugate(x) can have a non-zero imaginary part. And even if there is no noise, the type would be wrong: complex instead of real.

This causes the following doctest failure:

sage -t --long --warn-long 62.3 src/sage/symbolic/expression.pyx
**********************************************************************
File "src/sage/symbolic/expression.pyx", line 11105, in sage.symbolic.expression.Expression._plot_fast_callable
Failed example:
    plot(s)
Expected:
    Graphics object consisting of 1 graphics primitive
Got:
    Graphics object consisting of 0 graphics primitives
**********************************************************************

See pynac/pynac#117

Upstream: Reported upstream. Developers acknowledge bug.

CC: @rwst

Component: packages: standard

Reviewer: Jeroen Demeyer

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

[jdemeyer]

Branch: u/jdemeyer/ticket/19797

[vbraun]

comment:2

Isn't this sacrificing a lot of numerical performance? At least IEEE 754-2008 contains fused multiply-add.

The doctest is

            sage: x = var('x', domain='real')
            sage: s = abs((1+I*x)^4); s
            (I*x + 1)^2*(-I*x + 1)^2
            sage: s._plot_fast_callable(x)
            <sage.ext.interpreters.wrapper_py.Wrapper_py object at ...>
            sage: s._plot_fast_callable(x)(10)
            10201
            sage: abs((I*10+1)^4)
            10201
            sage: plot(s)
            Graphics object consisting of 1 graphics primitive

Isn't the real problem that abs does not default to hold=True

sage: ((1+I*x)^4).abs()
(I*x + 1)^2*(-I*x + 1)^2

so something that ought to be manifestly real is not. Compare:

sage: abs((1+I*x)^4)._plot_fast_callable(x)(RDF(10))
10201.00000000001
sage: type(_)
<type 'sage.rings.complex_double.ComplexDoubleElement'>

vs.

sage: ((1+I*x)^4).abs(hold=True)._plot_fast_callable(x)(RDF(10))
10201.000000000005
sage: type(_)
<type 'sage.rings.real_double.RealDoubleElement'>

---

New commits:

e468673	Simplify build of interpreters by skipping header files
694feb2	Build GSL in IEEE 754 compliant mode

[vbraun]

Commit: 694feb2

[jdemeyer]

comment:3

Replying to @vbraun:

Isn't this sacrificing a lot of numerical performance?

I have no idea but I doubt it. For example, Intel architectures added this instruction only very recently.

[vbraun]

comment:4

Yes but CPU's haven't seen any general speedups recently, either. Only more specialized instructions, and if you don't use them then you are not participating in progress.

[jdemeyer]

comment:5

Replying to @vbraun:

Yes but CPU's haven't seen any general speedups recently, either. Only more specialized instructions, and if you don't use them then you are not participating in progress.

I would say that at least in the specific case of complex multiplication, the use of FMA instructions is wrong.

[jdemeyer]

comment:6

Alternatively, we could probably play some tricks with volatile to fix complex multiplication.

[vbraun]

comment:7

So the question is: Should we guarantee that x*x.conj() has no numerical noise as imaginary part.

I'm not sure.

• Even besides FMA there is always the possibility of 80-bit x87 instructions messing things up, though thats being retired now.
• Blas are more than happy to use FMA, so if we make guarantees for scalars then they'll most likely be broken by diagonal matrices, say.

What I am pretty sure about is: multiplying x*x.conj() as complex numbers is a bad way of evaluating the absolute value. Fast callables using abs should return real and not complex values. And for that we need to hold abs. At which point the issue of the numerical noise in the imaginary part is pretty meaningless.

[rwst]

comment:8

Replying to @vbraun:

            sage: x = var('x', domain='real')
            sage: s = abs((1+I*x)^4); s
            (I*x + 1)^2*(-I*x + 1)^2

Isn't the real problem that abs does not default to hold=True

Maybe it's better to expand automatically if the result contains I. Automatic expansion also resolved a similar problem in #18952

sage: x = var('x', domain='real')
sage: ((1+I*x)^4).abs().expand()
x^4 + 2*x^2 + 1

so something that ought to be manifestly real is not.

It is but only when expanded.

[rwst]

comment:9

Maybe it's better to expand automatically if the result contains I.

Also only if the result is not of form abs(...)

[vbraun]

comment:10

And not of the form 1+abs(...) and so on

Automatically expanding isn't good enough either in general:

sage: x = var('x', domain='real')
sage: ((exp(I*x)+exp(-I*x))^4).abs().expand()
e^(4*I*x) + 4*e^(2*I*x) + 4*e^(-2*I*x) + e^(-4*I*x) + 6

[rwst]

Upstream: Reported upstream. Developers acknowledge bug.

[jdemeyer]

comment:12

Replying to @vbraun:

• Even besides FMA there is always the possibility of 80-bit x87 instructions messing things up

With x87, it is reasonable to assume that the whole imaginary part will be evaluated with 80 bits, so there won't be numerical noise either (see also the successful i386 buildbot reports).

[rwst]

comment:13

Replying to @vbraun:

And not of the form 1+abs(...) and so on

This can be fixed in GiNaC::abs_eval()

Automatically expanding isn't good enough either in general:

sage: x = var('x', domain='real')
sage: ((exp(I*x)+exp(-I*x))^4).abs().expand()
e^(4*I*x) + 4*e^(2*I*x) + 4*e^(-2*I*x) + e^(-4*I*x) + 6

Why not? This is real as well and could be rewritten as a sum of trig function expressions.

[jdemeyer]

comment:14

Replying to @rwst:

This is real as well and could be rewritten as a sum of trig function expressions.

Volker's point is about being "manifestly" real, which means: an expression which will always give a real result when evaluated, even if sub-expressions introduce numerical noise.

An expression like

e^(4*I*x) + e^(-4*I*x)

is not manifestly real, since there is no guarantee that the imaginary parts of the two terms will exactly cancel out.

[rwst]

comment:15

Well then, exclude functions in the treewalk in the search for I. Treewalks are very fast in comparison to other Pynac operations, and those are only a small fraction when compared to what Python costs.

[vbraun]

comment:16

Its not just that there is no guarantee against numerical noise in the imaginary part of exp(I*x)+exp(-I*x), the return type of the fast callable also changes from real to complex (comment:2). The symbolic expression containing abs() should only do automatic simplifications that are again manifestly real.

• Expanding polynomials over QQ[I] is safe, and that is the only safe case I can think of

• Transcendental functions are not safe

• Expanding polynomials whose coefficients are radical expressions over QQ[I] are not safe

Thoughts?

[rwst]

comment:17

Replying to @vbraun:

• Expanding polynomials whose coefficients are radical expressions over QQ[I] are not safe

Is this decidable without expansion? Simply weeding out rational exponents will make (I<sup>(5/4))*(I</sup>(3/4)) a false negative. OTOH first expanding then checking would be costly.

If so, I no longer think it useful.

[rwst]

comment:18

So, you're proposing to always hold abs() and probably have a keyword for x*conjugate(x) output?

[vbraun]

comment:19

Even with expansion I think it can be difficult to tell if a radical is real or not; Cardano's formula which can't be written in real radicals even if the roots are real comes to mind.

Of course it is possible to decide if a radical is real or not, for example using QQbar.

Can't we just abs(hold=True) by default? You can always simplify or set hold=False if you don't want that.

[jdemeyer]

comment:20

Replying to @vbraun:

Can't we just abs(hold=True) by default?

We could certainly do that as a stopgap-style solution. However, this is certainly going to break some doctests...

[rwst]

comment:21

Replying to @jdemeyer:

Replying to @vbraun:

Can't we just abs(hold=True) by default?

We could certainly do that as a stopgap-style solution. However, this is certainly going to break some doctests...

I get only two fails in functions/, symbolic/, calculus/, doc/ if I outcomment this snippet in Pynac which is responsible for the specific behaviour:
https://github.com/pynac/pynac/blob/master/ginac/inifcns.cpp#L253

[rwst]

comment:22

So the result is arrived at by two code parts: first abs of real power (or power of positive) is rewritten as power of abs and, second, even power of abs(x) is rewritten as x*x.conj() to the half power.

[rwst]

comment:23

So I'll just disable the abs(ex<sup>real)-->abs(ex)</sup>real part in that Pynac ticket if there are no objections.

[jdemeyer]

comment:24

Replying to @rwst:

second, even power of abs(x) is rewritten as x*x.conj() to the half power.

It's really this second part which is problematic.

[jdemeyer]

Changed dependencies from #19796 to none

[jdemeyer]

comment:27

The simplification abs(ex<sup>real)-->abs(ex)</sup>real is safe since the result is still manifestly real.

[rwst]

comment:28

This was introduced with pynac-0.3.9.2 from GiNaC. I'll revert it for pynac-0.5.4.

[jdemeyer]

Reviewer: Jeroen Demeyer

[jdemeyer]

comment:29

"Duplicate" of #19819.

[jdemeyer]

Changed author from Jeroen Demeyer to none

[jdemeyer]

Changed branch from u/jdemeyer/ticket/19797 to none

[jdemeyer]

Changed commit from 694feb2 to none