tight complex interval inverse #19964

miguelmarco · 2016-01-26T15:47:04Z

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

CC: @tscrim @videlec @dkrenn @jdemeyer @sagetrac-tmonteil @cheuberg

Component: numerical

Keywords: interval, root

Author: Miguel Marco

Branch: 39558b6

Reviewer: Vincent Delecroix, Marc Mezzarobba

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

The text was updated successfully, but these errors were encountered:

miguelmarco · 2016-01-26T15:48:12Z

Branch: u/mmarco/tight_complex_interval_inverse

sagetrac-git · 2016-01-26T15:55:59Z

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

`84ab655`	`compute tightly the inverse of complex interval.`
`a7a64ac`	`adapted division to new inverse`
`5b7082d`	`fixed doctests`
`97f4504`	`Merge branch 'inversecomplexintervals' into t/19964/tight_complex_interval_inverse`

sagetrac-git · 2016-01-26T15:55:59Z

Commit: 97f4504

miguelmarco · 2016-01-26T15:59:27Z

Author: Miguel Marco

miguelmarco · 2016-01-26T15:59:27Z

New commits:

`84ab655`	`compute tightly the inverse of complex interval.`
`a7a64ac`	`adapted division to new inverse`
`5b7082d`	`fixed doctests`
`97f4504`	`Merge branch 'inversecomplexintervals' into t/19964/tight_complex_interval_inverse`

miguelmarco · 2016-01-26T15:59:27Z

Changed keywords from none to interval, root

tscrim · 2016-01-26T16:21:45Z

comment:5

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

miguelmarco · 2016-01-26T17:31:51Z

comment:6

Added some examples to the description.

videlec · 2016-03-09T19:02:09Z

comment:7

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,

miguelmarco · 2016-03-09T19:24:12Z

comment:8

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

videlec · 2016-03-09T19:28:02Z

comment:9

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

miguelmarco · 2016-03-09T19:35:26Z

comment:10

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.

mezzarobba · 2016-03-25T10:41:55Z

comment:12

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

videlec · 2016-04-05T19:28:45Z

comment:13

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:

`1e5fe26`	`Trac 19964: add doctest for all branches`
`d57c626`	`Trac 19964: get rid of trailing whitespaces`

videlec · 2016-04-05T19:28:45Z

Reviewer: Vincent Delecroix

videlec · 2016-04-05T19:28:45Z

Changed commit from 97f4504 to d57c626

videlec · 2016-04-05T19:28:45Z

Changed branch from u/mmarco/tight_complex_interval_inverse to u/vdelecroix/19964

sagetrac-git · 2016-04-05T19:30:19Z

Changed commit from d57c626 to d57fee2

sagetrac-git · 2016-04-05T19:30:19Z

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

`421b3b3`	`Trac 19964: add doctest for all branches`
`d57fee2`	`Trac 19964: get rid of trailing whitespaces`

mezzarobba · 2016-04-06T08:10:21Z

comment:15

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

sagetrac-git · 2016-04-07T20:14:47Z

Changed commit from d57fee2 to 39558b6

sagetrac-git · 2016-04-07T20:14:47Z

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

`39558b6`	`Trac 19964: extend doctest`

mezzarobba · 2016-04-12T16:01:10Z

comment:17

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

mezzarobba · 2016-04-12T16:01:10Z

Changed reviewer from Vincent Delecroix to Vincent Delecroix, Marc Mezzarobba

vbraun · 2016-04-13T17:05:18Z

Changed branch from u/vdelecroix/19964 to 39558b6

jdemeyer · 2017-04-20T08:51:17Z

Changed commit from 39558b6 to none

jdemeyer · 2017-04-20T08:51:17Z

comment:19

Breakage: #22834

The tight reversion of complex intervals introduced in 84ab655 from sagemath#19964 made complex intervals wrong... and hence unreliable. This remained unnoticed until an apparently unrelated bug involving QQbar sagemath#37927. As an emergency measure we simply revert the tight reversion. Fixes sagemath#37927 ### 📝 Checklist  - [x] The title is concise and informative. - [x] The description explains in detail what this PR is about. - [x] I have linked a relevant issue or discussion. - [x] I have created tests covering the changes. - [x] I have updated the documentation and checked the documentation preview. URL: sagemath#37941 Reported by: Vincent Delecroix Reviewer(s): Marc Mezzarobba

…_ - thus for complex interval division. The implementation is inspired by Rokne-Lancaster 1971. Some history: - The first implementation of returning tight intervals was commit 84ba655 (see ticket sagemath#19964). - That implementation had bugs (returning too small intervals) and was reverted to fix ticket sagemath#37941. - This change is the second implementation.

This provides a (hopefully) correct implementation of computing a tight enclosure for the inverse of a complex interval. It fixes sagemath#37963. There is some history to this: Commit 84ab655 from sagemath#19964 introduced a tight enclosure but had some bugs such as ticket sagemath#37927. Thus, the that commit (and the subsequent partial fixes) were reverted in commit 8d59b12, see pull request sagemath#37941. This is a new implementation of Rokne-Lancaster that I wrote from scratch. ### 📝 Checklist - [x] The title is concise and informative. - [x] The description explains in detail what this PR is about. - [x] I have linked a relevant issue or discussion. - [x] I have created tests covering the changes. - [x] I have updated the documentation and checked the documentation preview. ### ⌛ Dependencies sagemath#37941 URL: sagemath#38360 Reported by: Matthias Goerner Reviewer(s): Marc Culler, Nathan Dunfield

miguelmarco added this to the sage-7.1 milestone Jan 26, 2016

miguelmarco added the p: major / 3 label Jan 26, 2016

This comment has been minimized.

Sign in to view

miguelmarco added s: needs review labels Jan 26, 2016

This comment has been minimized.

Sign in to view

videlec modified the milestones: sage-7.1, sage-7.2 Apr 5, 2016

mezzarobba added s: positive review and removed s: needs review labels Apr 12, 2016

vbraun removed the s: positive review label Apr 13, 2016

vbraun closed this as completed in db5a87e Apr 13, 2016

This was referenced May 5, 2024

right_kernel_matrix is wrong over QQbar #37927

Closed

Revert cif fast inversion 84ab655 #37941

Merged

unhyperbolic mentioned this issue Jul 14, 2024

Reintroducing tight complex interval inverse #38360

Merged

5 tasks

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

tight complex interval inverse #19964

tight complex interval inverse #19964

miguelmarco commented Jan 26, 2016

miguelmarco commented Jan 26, 2016

sagetrac-git mannequin commented Jan 26, 2016

sagetrac-git mannequin commented Jan 26, 2016

miguelmarco commented Jan 26, 2016

miguelmarco commented Jan 26, 2016

miguelmarco commented Jan 26, 2016

This comment has been minimized.

tscrim commented Jan 26, 2016

This comment has been minimized.

miguelmarco commented Jan 26, 2016

videlec commented Mar 9, 2016

miguelmarco commented Mar 9, 2016

videlec commented Mar 9, 2016

miguelmarco commented Mar 9, 2016

mezzarobba commented Mar 25, 2016

videlec commented Apr 5, 2016

videlec commented Apr 5, 2016

videlec commented Apr 5, 2016

videlec commented Apr 5, 2016

sagetrac-git mannequin commented Apr 5, 2016

sagetrac-git mannequin commented Apr 5, 2016

mezzarobba commented Apr 6, 2016

sagetrac-git mannequin commented Apr 7, 2016

sagetrac-git mannequin commented Apr 7, 2016

mezzarobba commented Apr 12, 2016

mezzarobba commented Apr 12, 2016

vbraun commented Apr 13, 2016

jdemeyer commented Apr 20, 2017

jdemeyer commented Apr 20, 2017

tight complex interval inverse #19964

tight complex interval inverse #19964

Comments

miguelmarco commented Jan 26, 2016

miguelmarco commented Jan 26, 2016

sagetrac-git mannequin commented Jan 26, 2016

sagetrac-git mannequin commented Jan 26, 2016

miguelmarco commented Jan 26, 2016

miguelmarco commented Jan 26, 2016

miguelmarco commented Jan 26, 2016

This comment has been minimized.

tscrim commented Jan 26, 2016

This comment has been minimized.

miguelmarco commented Jan 26, 2016

videlec commented Mar 9, 2016

miguelmarco commented Mar 9, 2016

videlec commented Mar 9, 2016

miguelmarco commented Mar 9, 2016

mezzarobba commented Mar 25, 2016

videlec commented Apr 5, 2016

videlec commented Apr 5, 2016

videlec commented Apr 5, 2016

videlec commented Apr 5, 2016

sagetrac-git mannequin commented Apr 5, 2016

sagetrac-git mannequin commented Apr 5, 2016

mezzarobba commented Apr 6, 2016

sagetrac-git mannequin commented Apr 7, 2016

sagetrac-git mannequin commented Apr 7, 2016

mezzarobba commented Apr 12, 2016

mezzarobba commented Apr 12, 2016

vbraun commented Apr 13, 2016

jdemeyer commented Apr 20, 2017

jdemeyer commented Apr 20, 2017