Verified computations

When used inside Sage, SnapPy can use interval arithmetic to rigorously verify the hyperbolicity of an orientable 3-manifold:

sage: M = Manifold("m015")
sage: M.verify_hyperbolicity()
(True,
 [0.6623589786224? + 0.5622795120623?*I,
  0.6623589786224? + 0.5622795120623?*I,
  0.6623589786224? + 0.5622795120623?*I])

It does so by producing intervals that are guaranteed to contain a true solution to the rectangular gluing equations (give bits_prec or dec_prec for higher precision intervals). This method works even when not all tetrahedra are positively oriented:

sage: M = Manifold("m015(3, 1)")
sage: M.tetrahedra_shapes('rect', intervals=True)
[0.625222762246? + 3.177940133813?*I,
 -0.0075523593782? + 0.5131157955971?*I,
 0.6515818912107? - 0.1955023488930?*I]

This is all achieved through a reimplementation of HIKMOT which pioneered the use of interval methods for hyperbolic manifolds. It can be used in a very similar way than HIKMOT, but uses Sage’s complex interval types and the Newton interval method (instead of the Krawczyk test) for certification. See Zgliczynski’s notes for a quick overview of these two tests.

This verification code was contributed by Matthias Goerner.

Methods of verifying hyperbolicity

snappy.verify.verify_hyperbolicity(manifold, verbose=False, bits_prec=53)

Given an orientable SnapPy Manifold, verifies its hyperbolicity. Similar to HIKMOT’s verify_hyperbolicity, the result is either (True, listOfShapeIntervals) or (False, []) if verification failed. listOfShapesIntervals is a list of complex intervals (elements in sage’s ComplexIntervalField) certified to contain the true shapes for the hyperbolic manifold.

Higher precision intervals can be obtained by setting bits_prec:

sage: from snappy import Manifold
sage: M = Manifold("m019")
sage: verify_hyperbolicity(M)
(True, [0.7805525278508? + 0.9144736629678?*I, 0.7805525278508? + 0.91447366296773?*I, 0.46002117557372? + 0.63262419360526?*I])

sage: M = Manifold("t02333(3,4)")
sage: verify_hyperbolicity(M)
(True, [2.152188153612? + 0.28494066790?*I, 1.92308491369? + 1.10360701507?*I, 0.014388591584? + 0.143084469681?*I, -2.5493670288? + 3.7453498408?*I, 0.142120333822? + 0.176540027036?*I, 0.50486686588? + 0.82829881681?*I, 0.50479249917? + 0.98036162786?*I, -0.58949570508? + 0.81267480427?*I])

A non-hyperbolic manifold:

sage: M = Manifold("t02333(1,0)")
sage: verify_hyperbolicity(M)
(False, [])

Under the hood, the function will call the CertifiedShapesEngine to produce intervals certified to contain a solution to the rectangular gluing equations. It then calls check_logarithmic_gluing_equations_and_positively_oriented_tets to verify that the logarithmic gluing equations are fulfilled and that all tetrahedra are positively oriented.

Generating certified shape intervals

The recommeded way to obtain certified intervals for the shapes is via manifold.tetrahedra_shapes(intervals=True) as described earlier. Here we document the CertifiedShapesEngine used internally to generate these intervals. It is of interest for those users who want to understand the underlying interval math and experiment with the Newton interval method.

class snappy.verify.CertifiedShapesEngine(M, initial_shapes, bits_prec=None, dec_prec=None)

An engine that is initialized with an approximated candidate solution to the rectangular gluing equations and produces intervals certified to contain a true solution. After the engine is successfully run, the resulting intervals are stored in certified_shapes which is a vector of elements in a Sage’s ComplexIntervalField.

A simple example to obtain certified shape intervals that uses CertifiedShapesEngine under the hood:

sage: from snappy import Manifold
sage: M = Manifold("m015")
sage: M.tetrahedra_shapes(bits_prec = 80, intervals = True)
[{'accuracies': (None, None, None, None), 'log': -0.140599787161480923256? + 0.703857721301476517492?*I, 'rect': 0.662358978622373012981? + 0.562279512062301243900?*I},
 {'accuracies': (None, None, None, None), 'log': -0.140599787161480923256? + 0.703857721301476517492?*I, 'rect': 0.662358978622373012981? + 0.562279512062301243900?*I},
 {'accuracies': (None, None, None, None), 'log': -0.140599787161480923256? + 0.703857721301476517492?*I, 'rect': 0.662358978622373012981? + 0.562279512062301243900?*I}]

Its objective is thus the same as HIKMOT and it is certainly HIKMOT inspired. However, it conceptually differs in that:

  1. It uses the Newton interval method instead of the Krawczyk test (we implement Gaussian elimination in interval arithmetic to compute the inverse of an interval matrix having interval arithmetic semantics, see mat_solve).

  2. It uses complex numbers in it’s Newton interval method. We simply use Sage’s complex interval type avoiding the need of converting n x n complex matrices into 2n x 2n real matrices as described Section 3.4 of the HIKMOT paper.

  3. We avoid automatic differentiation. We pick an independent set of equations of the following form and try to solve them:

    log(LHS) = 0

    where

    LHS = c * z0^a0 * (1-z0)^b0 * z1^a1 * (1-z1)^b1 * ...

    with a, b and c’s as returned by Manifold.gluing_equations(‘rect’).

    The derivative of log (LHS) with respect to zj is simply given by

    aj/zj - bj/(1-zj)

    and thus no need for automatic differentiation.

In contrast to HIKMOT, we use and return Sage’s native implementation of (complex) interval arithmetic here, which allows for increased interoperability. Another advantage is that Sage supports arbitrary precision. Unfortunately, performance suffers and this implementation is 5-10 times slower than HIKMOT.

Here is an example how to explicitly invoke the CertifiedShapesEngine:

sage: shapes = M.tetrahedra_shapes('rect', bits_prec = 80)
sage: C = CertifiedShapesEngine(M, shapes, bits_prec = 80)
sage: C.expand_until_certified()
True
sage: C.certified_shapes
(0.662358978622373012981? + 0.562279512062301243900?*I, 0.662358978622373012981? + 0.562279512062301243900?*I, 0.662358978622373012981? + 0.562279512062301243900?*I)
static certified_newton_iteration(equations, shapes)

Given shape intervals z, performs a Newton interval iteration N(z) as described in newton_iteration. Returns a pair (boolean, N(z)) where the boolean is True if N(z) is contained in z.

If the boolean is True, it is certified that N(z) contains a true solution, e.g., a point for which f is truely zero.

This follows from Theorem 1 of Zgliczynski’s notes.

Some examples:

sage: from snappy import Manifold
sage: M = Manifold("m019")
sage: C = CertifiedShapesEngine(M, M.tetrahedra_shapes('rect'),
...                           bits_prec = 80)

Intervals containing the true solution:

sage: good_shapes = vector([
...       C.CIF(C.RIF(0.78055, 0.78056), C.RIF(0.91447, 0.91448)),
...       C.CIF(C.RIF(0.78055, 0.78056), C.RIF(0.91447, 0.91448)),
...       C.CIF(C.RIF(0.46002, 0.46003), C.RIF(0.63262, 0.63263))])
sage: is_certified, shapes = CertifiedShapesEngine.certified_newton_iteration(C.equations, good_shapes)

sage: is_certified
True
sage: shapes
(0.78055253? + 0.91447366?*I, 0.78055253? + 0.91447367?*I, 0.46002118? + 0.63262420?*I)

This means that a true solution to the rectangular gluing equations is contained in both the given intervals (good_shapes) and the returned intervals (shapes) which are a refinement of the given intervals.

Intervals not containing a true solution:

sage: bad_shapes = vector([
...       C.CIF(C.RIF(0.78054, 0.78055), C.RIF(0.91447, 0.91448)),
...       C.CIF(C.RIF(0.78055, 0.78056), C.RIF(0.91447, 0.91448)),
...       C.CIF(C.RIF(0.46002, 0.46003), C.RIF(0.63262, 0.63263))])
sage: is_certified, shapes = CertifiedShapesEngine.certified_newton_iteration(C.equations, bad_shapes)
sage: is_certified
False
expand_until_certified(verbose=False)

Try to Newton interval iterate, then try to expand the shape intervals until we can certify they contain a true solution. If succeeded, return True and write certified shapes to certified_shapes. Set verbose = True for printing additional information.

static interval_vector_is_contained_in(vecA, vecB)

Given two vectors of intervals, return whether the first one is contained in the second one. Examples:

sage: RIF = RealIntervalField(80)
sage: CIF = ComplexIntervalField(80)
sage: box = CIF(RIF(-1,1),RIF(-1,1))
sage: a = [ CIF(0.1), CIF(1) + box ]
sage: b = [ CIF(0) + box, CIF(1) + 2 * box ]
sage: c = [ CIF(0), CIF(1) + 3 * box ]

sage: CertifiedShapesEngine.interval_vector_is_contained_in(a, b)
True
sage: CertifiedShapesEngine.interval_vector_is_contained_in(a, c)
False
sage: CertifiedShapesEngine.interval_vector_is_contained_in(b, a)
False
sage: CertifiedShapesEngine.interval_vector_is_contained_in(b, c)
False
sage: CertifiedShapesEngine.interval_vector_is_contained_in(c, a)
False
sage: CertifiedShapesEngine.interval_vector_is_contained_in(c, b)
False
static largest_diameter(shapes)

Given a vector of complex intervals, return the maximum of all their diameters:

sage: RIF = RealIntervalField(80)
sage: CIF = ComplexIntervalField(80)
sage: box = CIF(RIF(1,1.01),RIF(1,1.02))
sage: v = [ CIF(2) + box, CIF(3) + 2 * box ]
sage: CertifiedShapesEngine.largest_diameter(v)
0.019801980198019819393754
static log_gluing_LHS_derivatives(equations, shapes)

Compute the Jacobian of the vector-valued function f described in the above log_gluing_LHSs:

sage: from snappy import Manifold
sage: M = Manifold("m019")
sage: equations = M.gluing_equations('rect')
sage: RIF = RealIntervalField(80)
sage: CIF = ComplexIntervalField(80)
sage: shape1 = CIF(RIF(0.78055,0.78056), RIF(0.9144, 0.9145))
sage: shape2 = CIF(RIF(0.46002,0.46003), RIF(0.6326, 0.6327))
sage: shapes = [shape1, shape1, shape2]
sage: CertifiedShapesEngine.log_gluing_LHS_derivatives(equations, shapes)
[  0.292? - 1.667?*I   0.292? - 1.667?*I   0.752? - 1.034?*I]
[-0.5400? + 0.633?*I -0.5400? + 0.633?*I   1.561? + 1.829?*I]
[ 0.2482? + 1.034?*I  0.2482? + 1.034?*I  -2.313? - 0.795?*I]
[ 0.5400? - 0.633?*I -0.5400? + 0.633?*I                   0]
[-0.4963? - 2.068?*I  1.0800? - 1.266?*I   0.752? - 1.034?*I]
static log_gluing_LHSs(equations, shapes)

Given the result of M.gluing_equations(‘rect’) or a subset of rows of it and shapes, return a vector of log(LHS) where

LHS = c * z0 ** a0 * (1-z0) ** b0 * z1 ** a1 * ...

Let f: C^n -> C^n denote the function which takes shapes and returns the vector of log(LHS).

The reason we take the logarithm of the rectangular gluing equations is because the logarithmic derivative is of a particular nice form:

sage: from snappy import Manifold
sage: M = Manifold("m019")
sage: equations = M.gluing_equations('rect')
sage: RIF = RealIntervalField(80)
sage: CIF = ComplexIntervalField(80)
sage: zero = CIF(0).center()
sage: shape1 = CIF(RIF(0.78055,0.78056), RIF(0.9144, 0.9145))
sage: shape2 = CIF(RIF(0.46002,0.46003), RIF(0.6326, 0.6327))

An interval solution containing the true solution. The log of each rectangular equation should be 0 for the true solution, hence the interval should contain zero:

sage: shapes = [shape1, shape1, shape2]
sage: LHSs = CertifiedShapesEngine.log_gluing_LHSs(equations, shapes)
sage: LHSs
(0.000? + 0.000?*I, 0.000? + 0.000?*I, 0.000? + 0.000?*I, 0.000? + 0.000?*I, 0.000? + 0.000?*I)
sage: zero in LHSs[0]
True

An interval not containing the true solution:

sage: shapes = [shape1, shape1, shape1]
sage: LHSs = CertifiedShapesEngine.log_gluing_LHSs(equations, shapes)
sage: LHSs
(0.430? - 0.078?*I, -0.25? + 0.942?*I, -0.19? - 0.87?*I, 0.000? + 0.000?*I, 0.430? - 0.078?*I)
sage: zero in LHSs[0]
False
static mat_solve(m, v)

Given a matrix m and a vector v of (complex) intervals, returns the vector a such that v = m * a preserving interval arithmetics: if m’ is a matrix with values in the intervals of m and v’ is a vector with values in the intervals of v, then the intervals of the result a returned by this method are guarenteed to contain the entries of m’^-1 * v’.

Sage already provides a method for inverting matrices. However, it has a flaw and fails inverting interval matrices even though the interval determinant is far from containing zero (it returns unusuable matrices with entries (-inf, inf).

Our implementation improves on this by swapping rows to avoid diagonal entries close to zero during Gaussian elimination.

Setup a complex interval for example:

sage: RIF = RealIntervalField(80)
sage: CIF = ComplexIntervalField(80)
sage: fuzzy_four = CIF(RIF(3.9999,4.0001),RIF(-0.0001,0.0001))

Construct a matrix/vector with complex interval coefficients. One entry is a complex interval with non-zero diameter:

sage: m = matrix(CIF,
...      [  [ fuzzy_four, 3, 2, 3],
...         [          2, 3, 6, 2],
...         [          2, 4, 1, 6],
...         [          3, 2,-5, 2]])
sage: v = vector(CIF, [fuzzy_four, 2, 0 ,1])

Now compute the solutions a to v = m * a:

sage: a = CertifiedShapesEngine.mat_solve(m, v)
sage: a
(1.58? + 0.000?*I, -1.24? + 0.000?*I, 0.346? + 0.0000?*I, 0.24? + 0.000?*I)
sage: m * a
(4.0? + 0.00?*I, 2.0? + 0.00?*I, 0.0? + 0.00?*I, 1.00? + 0.00?*I)

The product actually contains the vector v, we check entry wise:

sage: [s in t for s, t in zip(v, m * a)]
[True, True, True, True]
static newton_iteration(equations, shapes)

Perform a Newton interval method of iteration for the function f described in log_gluing_LHSs.

Let z denote the shape intervals. Let z_center be a point close to the center point of the shape intervals (in the implementation, z_center is an interval of again, of length zero).

The result returned will be

N(z) = z_center - ((Df)(z))^-1 f(z_center)

A very approximate solution:

sage: from snappy import Manifold
sage: M = Manifold("m019")
sage: shapes = [ 0.7+1j, 0.7+1j, 0.5+0.5j ]

Get the equations and initialize zero-length intervals from it:

sage: C = CertifiedShapesEngine(M, shapes, bits_prec = 80)
sage: C.initial_shapes
(0.69999999999999995559107902? + 1*I, 0.69999999999999995559107902? + 1*I, 0.50000000000000000000000000? + 0.50000000000000000000000000?*I)

Do several Newton interval operations to get a better solution:

sage: shape_intervals = C.initial_shapes
sage: for i in range(4):
...     shape_intervals = CertifiedShapesEngine.newton_iteration(C.equations, shape_intervals)
...     print shape_intervals
(0.786746831183814577704? + 0.9208680745160821379529?*I, 0.7867468311838145777038? + 0.9208680745160821379529?*I, 0.4598680582870980309347? + 0.6194087185583516731751?*I)
(0.780561025176326485948? + 0.914496211844675048270?*I, 0.780561025176326485948? + 0.914496211844675048270?*I, 0.4599773577869384936554? + 0.632519407186945386957?*I)
(0.780552531045316100498? + 0.9144736621585220345231?*I, 0.7805525310453161004973? + 0.9144736621585220345231?*I, 0.4600211671037324947004? + 0.6326241909236695020810?*I)
(0.780552527850724832568? + 0.9144736629677264403330?*I, 0.7805525278507248325678? + 0.9144736629677264403330?*I, 0.4600211755737178641204? + 0.6326241936052562241142?*I)

For comparison:

sage: M.tetrahedra_shapes('rect')
[0.780552527850725 + 0.914473662967727*I, 0.780552527850725 + 0.914473662967726*I, 0.460021175573718 + 0.632624193605256*I]

Start with a rather big interval, note that the Newton interval method is stable in the sense that the interval size decreases:

sage: box = C.CIF(C.RIF(-0.0001,0.0001),C.RIF(-0.0001,0.0001))
sage: shape_intervals = C.initial_shapes.apply_map(lambda shape: shape + box)
sage: shape_intervals
(0.700? + 1.000?*I, 0.700? + 1.000?*I, 0.500? + 0.500?*I)
sage: for i in range(7):
...     shape_intervals = CertifiedShapesEngine.newton_iteration(C.equations, shape_intervals)
...     print shape_intervals
(0.79? + 0.92?*I, 0.79? + 0.92?*I, 0.460? + 0.62?*I)
(0.78? + 0.92?*I, 0.78? + 0.92?*I, 0.46? + 0.64?*I)
(0.781? + 0.915?*I, 0.7806? + 0.9145?*I, 0.4601? + 0.6327?*I)
(0.7805526? + 0.9144737?*I, 0.7805526? + 0.9144737?*I, 0.4600212? + 0.6326242?*I)
(0.780552527850725? + 0.914473662967727?*I, 0.780552527850725? + 0.914473662967727?*I, 0.4600211755737179? + 0.6326241936052562?*I)
(0.780552527850724837987? + 0.9144736629677264559386?*I, 0.7805525278507248379869? + 0.9144736629677264559386?*I, 0.4600211755737178728919? + 0.6326241936052561716379?*I)
(0.780552527850724837987? + 0.9144736629677264559386?*I, 0.7805525278507248379869? + 0.9144736629677264559386?*I, 0.4600211755737178728919? + 0.6326241936052561716379?*I)