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: |
Description (last modified by )
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 **********************************************************************
Change History (31)
comment:1 Changed 5 years ago by
- Branch set to u/jdemeyer/ticket/19797
comment:2 follow-ups: ↓ 3 ↓ 8 Changed 5 years ago by
- Cc rws added
- Commit set to 694feb2d71ebf521cf5f32dba69e0406d16e9f11
comment:3 in reply to: ↑ 2 Changed 5 years ago by
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: ↓ 5 Changed 5 years ago by
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
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.
comment:6 Changed 5 years ago by
Alternatively, we could probably play some tricks with volatile
to fix complex multiplication.
comment:7 follow-up: ↓ 12 Changed 5 years ago by
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: ↓ 9 Changed 5 years ago by
Replying to vbraun:
sage: x = var('x', domain='real') sage: s = abs((1+I*x)^4); s (I*x + 1)^2*(-I*x + 1)^2Isn'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
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: ↓ 13 Changed 5 years ago by
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
- 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
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: ↓ 14 Changed 5 years ago by
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
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.
comment:15 Changed 5 years ago by
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: ↓ 17 Changed 5 years ago by
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
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
So, you're proposing to always hold abs()
and probably have a keyword for x*conjugate(x)
output?
comment:19 follow-up: ↓ 20 Changed 5 years ago by
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: ↓ 21 Changed 5 years ago by
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
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: ↓ 24 Changed 5 years ago by
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
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
- 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 asx*x.conj()
to the half power.
It's really this second part which is problematic.
comment:25 Changed 5 years ago by
- Dependencies #19796 deleted
comment:26 Changed 5 years ago by
- Description modified (diff)
comment:27 Changed 5 years ago by
- 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
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
- 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
- Status changed from needs_review to positive_review
comment:31 Changed 5 years ago by
- Resolution set to duplicate
- Status changed from positive_review to closed
Isn't this sacrificing a lot of numerical performance? At least IEEE 754-2008 contains fused multiply-add.
The doctest is
Isn't the real problem that abs does not default to
hold=True
so something that ought to be manifestly real is not. Compare:
vs.
New commits:
Simplify build of interpreters by skipping header files
Build GSL in IEEE 754 compliant mode