Opened 5 years ago

Closed 5 years ago

#19797 closed enhancement (duplicate)

Disallow automatic abs() simplifications resulting non-real expressions

Reported by: jdemeyer Owned by:
Priority: major Milestone: sage-duplicate/invalid/wontfix
Component: packages: standard Keywords:
Cc: rws Merged in:
Authors: Reviewers: Jeroen Demeyer
Report Upstream: Reported upstream. Developers acknowledge bug. Work issues:
Branch: Commit:
Dependencies: Stopgaps:

Status badges

Description (last modified by 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 https://github.com/pynac/pynac/issues/117

Change History (31)

comment:1 Changed 5 years ago by jdemeyer

  • Branch set to u/jdemeyer/ticket/19797

comment:2 follow-ups: Changed 5 years ago by vbraun

  • Cc rws added
  • Commit set to 694feb2d71ebf521cf5f32dba69e0406d16e9f11

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:

e468673Simplify build of interpreters by skipping header files
694feb2Build GSL in IEEE 754 compliant mode
Last edited 5 years ago by vbraun (previous) (diff)

comment:3 in reply to: ↑ 2 Changed 5 years ago by jdemeyer

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.

comment:4 follow-up: Changed 5 years ago by 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.

comment:5 in reply to: ↑ 4 Changed 5 years ago by jdemeyer

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.

Last edited 5 years ago by jdemeyer (previous) (diff)

comment:6 Changed 5 years ago by jdemeyer

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

comment:7 follow-up: Changed 5 years ago by vbraun

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.

comment:8 in reply to: ↑ 2 ; follow-up: Changed 5 years ago by rws

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.

comment:9 in reply to: ↑ 8 Changed 5 years ago by rws

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

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

comment:10 follow-up: Changed 5 years ago by vbraun

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

comment:11 Changed 5 years ago by rws

  • Description modified (diff)
  • Report Upstream changed from N/A to Reported upstream. Developers acknowledge bug.

comment:12 in reply to: ↑ 7 Changed 5 years ago by jdemeyer

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).

comment:13 in reply to: ↑ 10 ; follow-up: Changed 5 years ago by rws

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.

comment:14 in reply to: ↑ 13 Changed 5 years ago by jdemeyer

Replying to rws:

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.

Last edited 5 years ago by jdemeyer (previous) (diff)

comment:15 Changed 5 years ago by rws

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.

comment:16 follow-up: Changed 5 years ago by vbraun

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?

comment:17 in reply to: ↑ 16 Changed 5 years ago by rws

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^(5/4))*(I^(3/4)) a false negative. OTOH first expanding then checking would be costly.

If so, I no longer think it useful.

comment:18 Changed 5 years ago by rws

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

comment:19 follow-up: Changed 5 years ago by vbraun

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.

comment:20 in reply to: ↑ 19 ; follow-up: Changed 5 years ago by 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...

comment:21 in reply to: ↑ 20 Changed 5 years ago by rws

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

comment:22 follow-up: Changed 5 years ago by rws

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.

comment:23 Changed 5 years ago by rws

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

comment:24 in reply to: ↑ 22 Changed 5 years ago by jdemeyer

  • Description modified (diff)
  • Summary changed from Build GSL in IEEE 754 compliant mode to Disallow automatic abs() simplifications resulting non-real expressions

Replying to rws:

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.

comment:25 Changed 5 years ago by jdemeyer

  • Dependencies #19796 deleted

comment:26 Changed 5 years ago by jdemeyer

  • Description modified (diff)

comment:27 Changed 5 years ago by jdemeyer

  • Description modified (diff)

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

comment:28 Changed 5 years ago by rws

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

comment:29 Changed 5 years ago by jdemeyer

  • Authors Jeroen Demeyer deleted
  • Branch u/jdemeyer/ticket/19797 deleted
  • Commit 694feb2d71ebf521cf5f32dba69e0406d16e9f11 deleted
  • Milestone changed from sage-7.0 to sage-duplicate/invalid/wontfix
  • Reviewers set to Jeroen Demeyer
  • Status changed from new to needs_review

"Duplicate" of #19819.

comment:30 Changed 5 years ago by jdemeyer

  • Status changed from needs_review to positive_review

comment:31 Changed 5 years ago by vbraun

  • Resolution set to duplicate
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.