Skip to content

Commit

Permalink
Trac #19964: tight complex interval inverse
Browse files Browse the repository at this point in the history
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)]
}}}

URL: http://trac.sagemath.org/19964
Reported by: mmarco
Ticket author(s): Miguel Marco
Reviewer(s): Vincent Delecroix, Marc Mezzarobba
  • Loading branch information
Release Manager authored and vbraun committed Apr 12, 2016
2 parents 0aec256 + 39558b6 commit db5a87e
Show file tree
Hide file tree
Showing 4 changed files with 233 additions and 75 deletions.
2 changes: 1 addition & 1 deletion src/sage/matrix/matrix0.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -1766,7 +1766,7 @@ cdef class Matrix(sage.structure.element.Matrix):
sage: K = (A-e).kernel()
sage: P = K.basis_matrix()
sage: P.str()
'[ 1.000000000000000? + 0.?e-17*I -2.116651487479748? + 0.0255565807096352?*I -0.2585224251020429? + 0.288602340904754?*I -0.4847545623533090? - 1.871890760086142?*I]'
'[ 1.000000000000000? + 0.?e-17*I -2.116651487479748? + 0.0255565807096352?*I -0.2585224251020429? + 0.2886023409047535?*I -0.4847545623533090? - 1.871890760086142?*I]'
Use single-row delimiters where appropriate::
Expand Down
44 changes: 22 additions & 22 deletions src/sage/matrix/matrix2.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -9044,23 +9044,23 @@ explicitly setting the argument to `True` or `False` will avoid this message."""
... [-1, 1, -6, -6, 5]])
sage: Q, R = A.QR()
sage: Q
[ -0.4588314677411235? -0.1260506983326509? 0.3812120831224489? -0.394573711338418? -0.687440062597?]
[ -0.4588314677411235? 0.4726901187474409? -0.05198346588033394? 0.717294125164660? -0.220962877263?]
[ 0.2294157338705618? 0.6617661662464172? 0.6619227988762521? -0.180872093737548? 0.1964114464561?]
[ 0.6882472016116853? 0.1890760474989764? -0.2044682991293135? 0.096630296654307? -0.662888631790?]
[ -0.4588314677411235? -0.1260506983326509? 0.3812120831224489? -0.394573711338418? -0.6874400625964?]
[ -0.4588314677411235? 0.4726901187474409? -0.05198346588033394? 0.7172941251646595? -0.2209628772631?]
[ 0.2294157338705618? 0.6617661662464172? 0.6619227988762521? -0.1808720937375480? 0.1964114464561?]
[ 0.6882472016116853? 0.1890760474989764? -0.2044682991293135? 0.0966302966543065? -0.6628886317894?]
[ -0.2294157338705618? 0.5357154679137663? -0.609939332995919? -0.536422031427112? 0.0245514308070?]
sage: R
[ 4.358898943540674? -0.4588314677411235? 13.07669683062202? 6.194224814505168? 2.982404540317303?]
[ 0 1.670171752907625? 0.5987408170800917? -1.292019657909672? 6.207996892883057?]
[ 0 0 5.444401659866974? 5.468660610611130? -0.682716185228386?]
[ 0 0 0 1.027626039419836? -3.61930014968662?]
[ 0 0 0 0 0.02455143080702?]
[ 0 0 5.444401659866974? 5.468660610611130? -0.6827161852283857?]
[ 0 0 0 1.027626039419836? -3.619300149686620?]
[ 0 0 0 0 0.024551430807012?]
sage: Q.conjugate_transpose()*Q
[1.000000000000000? 0.?e-18 0.?e-17 0.?e-15 0.?e-12]
[ 0.?e-18 1.000000000000000? 0.?e-16 0.?e-15 0.?e-12]
[ 0.?e-17 0.?e-16 1.000000000000000? 0.?e-15 0.?e-12]
[ 0.?e-15 0.?e-15 0.?e-15 1.000000000000000? 0.?e-12]
[ 0.?e-12 0.?e-12 0.?e-12 0.?e-12 1.000000000000?]
[1.000000000000000? 0.?e-18 0.?e-17 0.?e-16 0.?e-13]
[ 0.?e-18 1.000000000000000? 0.?e-17 0.?e-16 0.?e-13]
[ 0.?e-17 0.?e-17 1.000000000000000? 0.?e-16 0.?e-13]
[ 0.?e-16 0.?e-16 0.?e-16 1.000000000000000? 0.?e-13]
[ 0.?e-13 0.?e-13 0.?e-13 0.?e-13 1.0000000000000?]
sage: Q*R == A
True

Expand All @@ -9075,24 +9075,24 @@ explicitly setting the argument to `True` or `False` will avoid this message."""
sage: Q, R = A.QR()
sage: Q
[ -0.7302967433402215? 0.2070566455055649? + 0.5383472783144687?*I 0.2463049809998642? - 0.0764456358723292?*I 0.2381617683194332? - 0.1036596032779695?*I]
[ 0.0912870929175277? -0.2070566455055649? - 0.3778783780476559?*I 0.3786559533863032? - 0.1952221495524667?*I 0.701244450214469? - 0.364371165098660?*I]
[ 0.6390096504226938? + 0.0912870929175277?*I 0.1708217325420910? + 0.6677576817554466?*I -0.03411475806452072? + 0.04090198741767143?*I 0.3140171085506763? - 0.0825191718705412?*I]
[ 0.0912870929175277? -0.2070566455055649? - 0.3778783780476559?*I 0.3786559533863033? - 0.1952221495524667?*I 0.701244450214469? - 0.3643711650986595?*I]
[ 0.6390096504226938? + 0.0912870929175277?*I 0.1708217325420910? + 0.6677576817554466?*I -0.03411475806452072? + 0.04090198741767143?*I 0.3140171085506764? - 0.0825191718705412?*I]
[ 0.1825741858350554? + 0.0912870929175277?*I -0.03623491296347385? + 0.0724698259269477?*I 0.8632284069415110? + 0.06322839976356195?*I -0.4499694867611521? - 0.0116119181208918?*I]
sage: R
[ 10.95445115010333? 0.?e-18 - 1.917028951268082?*I 5.385938482134133? - 2.190890230020665?*I -0.2738612787525831? - 2.190890230020665?*I]
[ 0 4.829596256417300? + 0.?e-17*I -0.869637911123373? - 5.864879483945125?*I 0.993871898426712? - 0.3054085521207082?*I]
[ 0 4.829596256417300? + 0.?e-18*I -0.869637911123373? - 5.864879483945125?*I 0.993871898426712? - 0.3054085521207082?*I]
[ 0 0 12.00160760935814? + 0.?e-16*I -0.2709533402297273? + 0.4420629644486323?*I]
[ 0 0 0 1.942963944258992? + 0.?e-16*I]
sage: Q.conjugate_transpose()*Q
[1.000000000000000? + 0.?e-19*I 0.?e-18 + 0.?e-17*I 0.?e-17 + 0.?e-17*I 0.?e-16 + 0.?e-16*I]
[ 0.?e-18 + 0.?e-17*I 1.000000000000000? + 0.?e-17*I 0.?e-17 + 0.?e-17*I 0.?e-16 + 0.?e-16*I]
[ 0.?e-17 + 0.?e-17*I 0.?e-17 + 0.?e-17*I 1.000000000000000? + 0.?e-16*I 0.?e-16 + 0.?e-16*I]
[ 0.?e-16 + 0.?e-16*I 0.?e-16 + 0.?e-16*I 0.?e-16 + 0.?e-16*I 1.000000000000000? + 0.?e-15*I]
[ 0.?e-17 + 0.?e-17*I 0.?e-17 + 0.?e-17*I 1.000000000000000? + 0.?e-17*I 0.?e-16 + 0.?e-16*I]
[ 0.?e-16 + 0.?e-16*I 0.?e-16 + 0.?e-16*I 0.?e-16 + 0.?e-16*I 1.000000000000000? + 0.?e-16*I]
sage: Q*R - A
[ 0.?e-17 0.?e-17 + 0.?e-17*I 0.?e-16 + 0.?e-16*I 0.?e-16 + 0.?e-16*I]
[ 0.?e-18 0.?e-17 + 0.?e-17*I 0.?e-16 + 0.?e-16*I 0.?e-15 + 0.?e-15*I]
[ 0.?e-18 0.?e-17 + 0.?e-17*I 0.?e-16 + 0.?e-16*I 0.?e-16 + 0.?e-16*I]
[0.?e-17 + 0.?e-18*I 0.?e-17 + 0.?e-17*I 0.?e-16 + 0.?e-16*I 0.?e-16 + 0.?e-16*I]
[0.?e-18 + 0.?e-18*I 0.?e-18 + 0.?e-17*I 0.?e-16 + 0.?e-16*I 0.?e-15 + 0.?e-16*I]
[0.?e-18 + 0.?e-18*I 0.?e-18 + 0.?e-18*I 0.?e-16 + 0.?e-16*I 0.?e-16 + 0.?e-16*I]

A rank-deficient rectangular matrix, with both values of the ``full`` keyword. ::

Expand Down Expand Up @@ -10495,9 +10495,9 @@ explicitly setting the argument to `True` or `False` will avoid this message."""
True
sage: _, T = A.is_similar(B, transformation=True)
sage: T
[ 1.0000000000000? + 0.?e-13*I 0.?e-13 + 0.?e-13*I 0.?e-13 + 0.?e-13*I]
[-0.6666666666667? + 0.?e-13*I 0.16666666666667? + 0.?e-14*I -0.8333333333334? + 0.?e-13*I]
[ 0.6666666666667? + 0.?e-13*I 0.?e-13 + 0.?e-13*I -0.333333333334? + 0.?e-13*I]
[ 1.00000000000000? + 0.?e-14*I 0.?e-14 + 0.?e-14*I 0.?e-14 + 0.?e-14*I]
[-0.66666666666667? + 0.?e-15*I 0.166666666666667? + 0.?e-15*I -0.83333333333334? + 0.?e-14*I]
[ 0.66666666666667? + 0.?e-14*I 0.?e-14 + 0.?e-14*I -0.33333333333333? + 0.?e-14*I]
sage: T.change_ring(QQ)
[ 1 0 0]
[-2/3 1/6 -5/6]
Expand Down
Loading

0 comments on commit db5a87e

Please sign in to comment.