Opened 5 years ago

Closed 5 years ago

Last modified 4 years ago

#19964 closed enhancement (fixed)

tight complex interval inverse

Reported by: mmarco Owned by:
Priority: major Milestone: sage-7.2
Component: numerical Keywords: interval, root
Cc: tscrim, vdelecroix, dkrenn, jdemeyer, tmonteil, cheuberg Merged in:
Authors: Miguel Marco Reviewers: Vincent Delecroix, Marc Mezzarobba
Report Upstream: N/A Work issues:
Branch: 39558b6 (Commits) Commit:
Dependencies: Stopgaps:

Description (last modified by mmarco)

The current method to compute the multiplicative inverse of a complex interval uses schoolbok formula that produces intervals much bigger than the actual result.

This ticket implements a method that produces a tight interval enclosure of the result.

As a result, we get better precission in root isolating methods (and hence, less steps needed).

Here are some examples that illustrate the impact of the changes:

Old behaviour:

sage: a = CIF((1, 2), (3, 4))
sage: a.real().lower(), a.real().upper()
(1.00000000000000, 2.00000000000000)
sage: a = CIF((1, 2), (3, 4))
sage: b = a.__invert__()
sage: b.real().lower(), b.real().upper()
(0.0499999999999999, 0.200000000000001)
sage: b.imag().lower(), b.imag().upper()
(-0.400000000000001, -0.149999999999999)
sage: b.diameter()
1.20000000000000
sage: %time a.__invert__()
CPU times: user 0 ns, sys: 0 ns, total: 0 ns
Wall time: 23.1 µs
1.? - 1.?*I

new behaviour:

sage: a = CIF((1, 2), (3, 4))
sage: a.real().lower(), a.real().upper()
(1.00000000000000, 2.00000000000000)
sage: b = a.__invert__()
sage: b.real().lower(), b.real().upper()
(0.0588235294117647, 0.153846153846154)
sage: b.imag().lower(), b.imag().upper()
(-0.300000000000001, -0.200000000000000)
sage: b.diameter()
0.893617021276596
sage: %time a.__invert__()
CPU times: user 0 ns, sys: 0 ns, total: 0 ns
Wall time: 19.1 µs
0.1? - 0.3?*I

Another example with a smaller interval,

old:

sage: a = CIF((5, 5.01), (7, 7.01))
sage: b = a.__invert__()
sage: b.diameter()
0.00523867991752868
sage: b.real().lower(), b.real().upper()
(0.0673489564952680, 0.0677027027027028)
sage: b.imag().lower(), b.imag().upper()
(-0.0947297297297298, -0.0942885390933752)
sage: %time a.__invert__()
CPU times: user 0 ns, sys: 0 ns, total: 0 ns
Wall time: 23.1 µs
-0.068? - 0.095?*I

new:

sage: a = CIF((5, 5.01), (7, 7.01))
sage: b = a.__invert__()
sage: b.diameter()
0.00253766599329326
sage: b.real().lower(), b.real().upper()
(0.0674398874563158, 0.0676112447891434)
sage: b.imag().lower(), b.imag().upper()
(-0.0945945945945946, -0.0944232370063658)
sage: a = CIF((-5.01, -5), (7, 7.01))
sage: %time a.__invert__()
CPU times: user 0 ns, sys: 0 ns, total: 0 ns
Wall time: 49.1 µs
-0.068? - 0.0945?*I

As you can see, the diameter of the result can easily get cut to half or even smaller.

The timings of the inversion can vary, but they remain in the same order of magnitude as before. It has a noticeable effect in root isolation:

old:

sage: R.<x> = QQ[]
sage: f = -7/6*x^7 - x^6 + 1/2*x^4 + 2/17*x^3 + 2*x^2 - 6*x - 9
sage: %time f.roots(CC)
CPU times: user 6 ms, sys: 0 ns, total: 6 ms
Wall time: 5.34 ms
[(-1.03157268729039, 1),
 (-1.18964736821330 - 0.850802715570186*I, 1),
 (-1.18964736821330 + 0.850802715570186*I, 1),
 (0.0780792982605128 - 1.39288589462175*I, 1),
 (0.0780792982605128 + 1.39288589462175*I, 1),
 (1.19878298502655 - 0.599304323571307*I, 1),
 (1.19878298502655 + 0.599304323571307*I, 1)]
sage: %time f.roots(QQbar)
CPU times: user 17 ms, sys: 0 ns, total: 17 ms
Wall time: 15.7 ms
[(-1.031572687290387?, 1),
 (-1.189647368213303? - 0.8508027155701860?*I, 1),
 (-1.189647368213303? + 0.8508027155701860?*I, 1),
 (0.07807929826051275? - 1.392885894621755?*I, 1),
 (0.07807929826051275? + 1.392885894621755?*I, 1),
 (1.198782985026555? - 0.5993043235713073?*I, 1),
 (1.198782985026555? + 0.5993043235713073?*I, 1)]

new:

sage: R.<x>=QQ[]
sage: f = -7/6*x^7 - x^6 + 1/2*x^4 + 2/17*x^3 + 2*x^2 - 6*x - 9
sage: %time f.roots(CC)
CPU times: user 4 ms, sys: 0 ns, total: 4 ms
Wall time: 3.28 ms
[(-1.03157268729039, 1),
 (-1.18964736821330 - 0.850802715570186*I, 1),
 (-1.18964736821330 + 0.850802715570186*I, 1),
 (0.0780792982605128 - 1.39288589462175*I, 1),
 (0.0780792982605128 + 1.39288589462175*I, 1),
 (1.19878298502655 - 0.599304323571307*I, 1),
 (1.19878298502655 + 0.599304323571307*I, 1)]
sage: %time f.roots(QQbar)
CPU times: user 9 ms, sys: 1 ms, total: 10 ms
Wall time: 9.3 ms
[(-1.031572687290387?, 1),
 (-1.189647368213303? - 0.8508027155701860?*I, 1),
 (-1.189647368213303? + 0.8508027155701860?*I, 1),
 (0.07807929826051275? - 1.392885894621755?*I, 1),
 (0.07807929826051275? + 1.392885894621755?*I, 1),
 (1.198782985026555? - 0.5993043235713073?*I, 1),
 (1.198782985026555? + 0.5993043235713073?*I, 1)]

Change History (19)

comment:1 Changed 5 years ago by mmarco

  • Branch set to u/mmarco/tight_complex_interval_inverse

comment:2 Changed 5 years ago by git

  • Commit set to 97f4504c20b5bfbf315f5e00013dab27df25fbd3

Branch pushed to git repo; I updated commit sha1. New commits:

84ab655compute tightly the inverse of complex interval.
a7a64acadapted division to new inverse
5b7082dfixed doctests
97f4504Merge branch 'inversecomplexintervals' into t/19964/tight_complex_interval_inverse

comment:3 Changed 5 years ago by mmarco

  • Authors set to Miguel Marco
  • Cc tscrim vdelecroix dkren jdemeyern added
  • Component changed from PLEASE CHANGE to numerical
  • Description modified (diff)
  • Keywords interval root added
  • Status changed from new to needs_review
  • Type changed from PLEASE CHANGE to enhancement

New commits:

84ab655compute tightly the inverse of complex interval.
a7a64acadapted division to new inverse
5b7082dfixed doctests
97f4504Merge branch 'inversecomplexintervals' into t/19964/tight_complex_interval_inverse

comment:4 Changed 5 years ago by tscrim

  • Cc jdemeyer added; jdemeyern removed

comment:5 Changed 5 years ago by tscrim

  • Cc dkrenn added; dkren removed

Can you give an example which clearly illustrates the better bounds? Also, have you run timings versus the schoolbook version? (I don't have the technical abilities to review this ticket though.)

comment:6 Changed 5 years ago by mmarco

  • Description modified (diff)

Added some examples to the description.

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

comment:7 Changed 5 years ago by vdelecroix

Is this better in all cases?

-             0.?e-182 - 0.7598602580415435?*I,
-             0.?e-249 + 0.7598602580415435?*I,
+             0.?e-215 - 0.7598602580415435?*I,
+             0.?e-229 + 0.7598602580415435?*I,

comment:8 Changed 5 years ago by mmarco

It should be, up to floating point artifacts. I don't know much about the subtleties of floating point numbers though.

comment:9 Changed 5 years ago by vdelecroix

I was wondering about the e-249 becoming a e-229 in the above example...

comment:10 Changed 5 years ago by mmarco

Yes, that is what i imagined, but i would guess that in those two cases the discrepancy is due to those pesky numerical artifacts that can appear just by, for example, performing the same operations in different order.

comment:11 Changed 5 years ago by tscrim

  • Cc tmonteil cheuberg added

comment:12 Changed 5 years ago by mmezzarobba

IMHO the tests should cover all branches of the new implementation (ideally with test cases where the rounding modes matter, if that's not too hard to achieve).

comment:13 Changed 5 years ago by vdelecroix

  • Branch changed from u/mmarco/tight_complex_interval_inverse to u/vdelecroix/19964
  • Commit changed from 97f4504c20b5bfbf315f5e00013dab27df25fbd3 to d57c62632a612d8f098515ff604ba153a3566ec2
  • Milestone changed from sage-7.1 to sage-7.2
  • Reviewers set to Vincent Delecroix

Following comment:12 I added a doctest. And I also removed the trailing whitespaces in a second commit.

This is good for me to go. Any remarks?


New commits:

1e5fe26Trac 19964: add doctest for all branches
d57c626Trac 19964: get rid of trailing whitespaces

comment:14 Changed 5 years ago by git

  • Commit changed from d57c62632a612d8f098515ff604ba153a3566ec2 to d57fee2259bf69d0f0471f1bc89c3a81f1e12e27

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

421b3b3Trac 19964: add doctest for all branches
d57fee2Trac 19964: get rid of trailing whitespaces

comment:15 Changed 5 years ago by mmezzarobba

lgtm except that I'd suggest the following change in order to really test all possible cases:

@@ -1130,12 +1130,15 @@ cdef class ComplexIntervalFieldElement(sage.structure.element.FieldElement):
         Check that the code is valid in all regions of the complex plane::
 
             sage: rats = [15423/906, 337/59976, 145151/145112]
-            sage: for a in rats:
-            ....:     for b in rats:
-            ....:         x = CIF(a, b)
-            ....:         assert (x * (~x) - 1).contains_zero()
-            ....:         x = -CIF(a,b)
-            ....:         assert (x * (~x) - 1).contains_zero()
+            sage: ivs = [RIF(aa, bb) for a in rats for b in rats
+            ....:                   for aa in [-a, a] for bb in [-b, b]
+            ....:                   if aa <= bb]
+            sage: for x in ivs:
+            ....:     for y in ivs:
+            ....:         z = CIF(x, y)
+            ....:         inv = ~z
+            ....:         assert ((z*inv - 1).contains_zero()
+            ....:                 or z.contains_zero() and inv.real().is_NaN())

comment:16 Changed 5 years ago by git

  • Commit changed from d57fee2259bf69d0f0471f1bc89c3a81f1e12e27 to 39558b6299664e10bd008e2f97bb9c64e62096b3

Branch pushed to git repo; I updated commit sha1. New commits:

39558b6Trac 19964: extend doctest

comment:17 Changed 5 years ago by mmezzarobba

  • Reviewers changed from Vincent Delecroix to Vincent Delecroix, Marc Mezzarobba
  • Status changed from needs_review to positive_review

You really don't want to test the cases where the interval contains zero ;-). But so be it.

comment:18 Changed 5 years ago by vbraun

  • Branch changed from u/vdelecroix/19964 to 39558b6299664e10bd008e2f97bb9c64e62096b3
  • Resolution set to fixed
  • Status changed from positive_review to closed

comment:19 Changed 4 years ago by jdemeyer

  • Commit 39558b6299664e10bd008e2f97bb9c64e62096b3 deleted

Breakage: #22834

Note: See TracTickets for help on using tickets.