docs missing on the question of precision

[dimpase]

We're trying to use the package for computing fundamental domains of a discrete subgroup of $SL_4(\mathbb{R})$, and wonder about the precision issues. There will be singularities (many more hypersurfaces, say, 10 or 15, intersecting in a point, than in the generic case, i.e. 4).
Can such a case be dealt with directly?

Does it use the machine precision? Or something better? The docs are silent on this.
@Ilia-Smilga

[PBrdng]

Thanks for reaching out. In principle, singularities should not pose a problem, since the algorithm is computing critical points of $\log\ (f_1\cdots f_k)$, where the $f_i$ define your hypersurfaces. Of course, there can be other numerical issues.

What is the number of hypersurfaces and how many variables are we talking about?

The implementation uses 64-bit machine precision by using methods from HomotopyContinuation.jl and DifferentialEquations.jl out-of-the-box. That is why we didn't document it.

Can you share a bit more details on your example, maybe I can also take a look at it.

[dimpase]

The background of the computations we'd like to carry out is geometric group theory, more precisely, computing a fundamental domain of the Hilbert modular group of a real quadratic extension $\mathbb{K}$ of $\mathbb{Q}$ - this is
$SL_2(O_{\mathbb{K}})$, where $O_\mathbb{K}$ is the ring of integers of $\mathbb{K}$. (And same for some of its subgroups). Such domains exist in $\mathbb{R}^4$ - specifically in the product of two copies of the hyperbolic plane.

The hypersurfaces are mostly affine of degree 4 (a particular rational type of them, although as far as I understand this is of no direct help for your solver). (there are also few of them of smaller degree - in particular e.g. $y_1$ and $y_2$ are always positive), with coefficients in $\mathbb{K}$ (we can get rid of the latter, by making formulae messier).
Here is a typical part of such input - here we have $\mathbb{K}=\mathbb{Q}[\sqrt{3}]$:

@var x1, x2, y1, y2
w=sqrt(3)
p=[w^4 - 4*w^3*x1 + 4*w^2*x1^2 + 4*w^3*x2 - 16*w^2*x1*x2 + 16*w*x1^2*x2 + 4*w^2*x2^2 - 16*w*x1*x2^2 + 16*x1^2*x2^2 + 4*w^2*y1^2 + 16*w*x2*y1^2 + 16*x2^2*y1^2 + 4*w^2*y2^2 - 16*w*x1*y2^2 + 16*x1^2*y2^2 + 16*y1^2*y2^2 - 1,
        w^4*x1^2*x2^2 + w^4*x2^2*y1^2 + w^4*x1^2*y2^2 + w^4*y1^2*y2^2 + 2*w^3*x1^2*x2 - 2*w^3*x1*x2^2 - 2*w^2*x1^2*x2^2 + 2*w^3*x2*y1^2 - 2*w^2*x2^2*y1^2 - 2*w^3*x1*y2^2 - 2*w^2*x1^2*y2^2 - 2*w^2*y1^2*y2^2 + 2*w^2*x1^2*x2 + 2*w^2*x1*x2^2 + 2*w^2*x2*y1^2 + 2*w^2*x1*y2^2 + w^2*x1^2 - 4*w^2*x1*x2 - 2*w*x1^2*x2 + w^2*x2^2 + 2*w*x1*x2^2 + x1^2*x2^2 + w^2*y1^2 - 2*w*x2*y1^2 + x2^2*y1^2 + w^2*y2^2 + 2*w*x1*y2^2 + x1^2*y2^2 + y1^2*y2^2 + 2*w*x1^2 - 2*x1^2*x2 - 2*w*x2^2 - 2*x1*x2^2 + 2*w*y1^2 - 2*x2*y1^2 - 2*w*y2^2 - 2*x1*y2^2 - 2*w*x1 + x1^2 + 2*w*x2 + 4*x1*x2 + x2^2 + y1^2 + y2^2 - 2*x1 - 2*x2,
        x1^2*x2^2 + x2^2*y1^2 + x1^2*y2^2 + y1^2*y2^2 + 2*x1^2*x2 + 2*x1*x2^2 + 2*x2*y1^2 + 2*x1*y2^2 + x1^2 + 4*x1*x2 + x2^2 + y1^2 + y2^2 + 2*x1 + 2*x2,
        w^4*x1^2*x2^2 + w^4*x2^2*y1^2 + w^4*x1^2*y2^2 + w^4*y1^2*y2^2 + 2*w^3*x1^2*x2 - 2*w^3*x1*x2^2 + 2*w^3*x2*y1^2 - 2*w^3*x1*y2^2 + w^2*x1^2 - 4*w^2*x1*x2 + w^2*x2^2 + w^2*y1^2 + w^2*y2^2 - 2*w*x1 + 2*w*x2,
        w^4 - 2*w^3*x1 + w^2*x1^2 + 2*w^3*x2 - 4*w^2*x1*x2 + 2*w*x1^2*x2 + w^2*x2^2 - 2*w*x1*x2^2 + x1^2*x2^2 + w^2*y1^2 + 2*w*x2*y1^2 + x2^2*y1^2 + w^2*y2^2 - 2*w*x1*y2^2 + x1^2*y2^2 + y1^2*y2^2 - 2*w^2*x1 + 2*w*x1^2 - 2*w^2*x2 + 2*x1^2*x2 - 2*w*x2^2 + 2*x1*x2^2 + 2*w*y1^2 + 2*x2*y1^2 - 2*w*y2^2 + 2*x1*y2^2 - 2*w^2 + 2*w*x1 + x1^2 - 2*w*x2 + 4*x1*x2 + x2^2 + y1^2 + y2^2 + 2*x1 + 2*x2,
        w^4*x1^2*x2^2 + w^4*x2^2*y1^2 + w^4*x1^2*y2^2 + w^4*y1^2*y2^2 + 4*w^3*x1^2*x2 - 4*w^3*x1*x2^2 + 4*w^3*x2*y1^2 - 4*w^3*x1*y2^2 + 4*w^2*x1^2 - 16*w^2*x1*x2 + 4*w^2*x2^2 + 4*w^2*y1^2 + 4*w^2*y2^2 - 16*w*x1 + 16*w*x2 + 15,
        x1^2*x2^2 + x2^2*y1^2 + x1^2*y2^2 + y1^2*y2^2 - 1,
        w^4 - 2*w^3*x1 + w^2*x1^2 + 2*w^3*x2 - 4*w^2*x1*x2 + 2*w*x1^2*x2 + w^2*x2^2 - 2*w*x1*x2^2 + x1^2*x2^2 + w^2*y1^2 + 2*w*x2*y1^2 + x2^2*y1^2 + w^2*y2^2 - 2*w*x1*y2^2 + x1^2*y2^2 + y1^2*y2^2 + 2*w^2*x1 - 2*w*x1^2 + 2*w^2*x2 - 2*x1^2*x2 + 2*w*x2^2 - 2*x1*x2^2 - 2*w*y1^2 - 2*x2*y1^2 + 2*w*y2^2 - 2*x1*y2^2 - 2*w^2 + 2*w*x1 + x1^2 - 2*w*x2 + 4*x1*x2 + x2^2 + y1^2 + y2^2 - 2*x1 - 2*x2,
        x1^2*x2^2 + x2^2*y1^2 + x1^2*y2^2 + y1^2*y2^2 - 2*x1^2*x2 - 2*x1*x2^2 - 2*x2*y1^2 - 2*x1*y2^2 + x1^2 + 4*x1*x2 + x2^2 + y1^2 + y2^2 - 2*x1 - 2*x2,
        w^4*x1^2*x2^2 + w^4*x2^2*y1^2 + w^4*x1^2*y2^2 + w^4*y1^2*y2^2 - 2*w^4*x1^2*x2 - 2*w^4*x1*x2^2 - 2*w^4*x2*y1^2 - 2*w^4*x1*y2^2 + w^4*x1^2 + 4*w^4*x1*x2 + 2*w^3*x1^2*x2 + w^4*x2^2 - 2*w^3*x1*x2^2 - 2*w^2*x1^2*x2^2 + w^4*y1^2 + 2*w^3*x2*y1^2 - 2*w^2*x2^2*y1^2 + w^4*y2^2 - 2*w^3*x1*y2^2 - 2*w^2*x1^2*y2^2 - 2*w^2*y1^2*y2^2 - 2*w^4*x1 - 2*w^3*x1^2 - 2*w^4*x2 + 6*w^2*x1^2*x2 + 2*w^3*x2^2 + 6*w^2*x1*x2^2 - 2*w^3*y1^2 + 6*w^2*x2*y1^2 + 2*w^3*y2^2 + 6*w^2*x1*y2^2 + w^4 + 2*w^3*x1 - 3*w^2*x1^2 - 2*w^3*x2 - 20*w^2*x1*x2 - 2*w*x1^2*x2 - 3*w^2*x2^2 + 2*w*x1*x2^2 + x1^2*x2^2 - 3*w^2*y1^2 - 2*w*x2*y1^2 + x2^2*y1^2 - 3*w^2*y2^2 + 2*w*x1*y2^2 + x1^2*y2^2 + y1^2*y2^2 + 12*w^2*x1 + 4*w*x1^2 + 12*w^2*x2 - 4*x1^2*x2 - 4*w*x2^2 - 4*x1*x2^2 + 4*w*y1^2 - 4*x2*y1^2 - 4*w*y2^2 - 4*x1*y2^2 - 8*w^2 - 8*w*x1 + 4*x1^2 + 8*w*x2 + 16*x1*x2 + 4*x2^2 + 4*y1^2 + 4*y2^2 - 16*x1 - 16*x2 + 15,
        w^4 - 2*w^3*x1 + w^2*x1^2 + 2*w^3*x2 - 4*w^2*x1*x2 + 2*w*x1^2*x2 + w^2*x2^2 - 2*w*x1*x2^2 + x1^2*x2^2 + w^2*y1^2 + 2*w*x2*y1^2 + x2^2*y1^2 + w^2*y2^2 - 2*w*x1*y2^2 + x1^2*y2^2 + y1^2*y2^2 - 1]

This particular fragment features a singularity - a common root at the centre of the ball of squared radius 3/5 inside the sphere given by the equation

(x1-sqrt(3)/2)^2+(x2+sqrt(3)/2)^2+(y1-0.5)^2+(y2-0.5)^2-0.6

The full input is typically much longer - however we're programming a subdivision-based algorithm allowing to "zoom in" on subsets with (hopefully) relatively few hypersurfaces (tens or hundreds).

We are only interested in the cell where all the functions (specified by the hypersurface equations) are positive, and the adjacent cells (where just one function is negative).
Does the solver allow such restricted problems, or it's not something that the underlying algorithm can deal with?

It is crucial to us, as we cannot lose any topological information due to numerical issues (the topological information is needed for computing a presentation of the group, which is the real problem we're after), does the solver detect such issues and throws an error? And related question - is there a way to get a correctness certificate for the computation (something that's often avaliable in convex optimisation)?

Last but not the least - I presume that an incremental (by adding more hyperplanes) refining of the arrangement isn't feasible with this approach, is it true?

[dimpase]

Does the solver allow such restricted problems

one can imagine that at least after the critical points are found, only the ones within the desired region are taken for further processing, and the rest ignored.

[PBrdng]

We are only interested in the cell where all the functions (specified by the hypersurface equations) are positive, and the adjacent cells (where just one function is negative).
Does the solver allow such restricted problems, or it's not something that the underlying algorithm can deal with?

The algorithm works in two stages. In the first stage, the critical points of $\log (f_1\cdots f_k)$ are found. In the second stage, the algorithm connects the critical points that lie in the same connected component. In principle, one can restrict the second stage so that it only processes cells where all functions are positive. However, the second stage is usually the stage which is computationally much less demanding. The bottleneck is really the first stage and I don't see a way to restrict it.

It is crucial to us, as we cannot lose any topological information due to numerical issues (the topological information is needed for computing a presentation of the group, which is the real problem we're after), does the solver detect such issues and throws an error?

There are a lot of built-in heuristics to detect numerical issues. I would trust the numerical computation, but of course it is not a certified computation.

And related question - is there a way to get a correctness certificate for the computation (something that's often avaliable in convex optimisation)?

The answer to this question is subtle. It is possible to certify the computation of each critical point individually (although, this is currently not implemented as an option here). However, there is no method (not even theoretically) to certify whether we have found all critical points.

Last but not the least - I presume that an incremental (by adding more hyperplanes) refining of the arrangement isn't feasible with this approach, is it true?

This is true, because adding more hyperplanes will change the critical points, so the previous computation is useless.

one can imagine that at least after the critical points are found, only the ones within the desired region are taken for further processing, and the rest ignored.

See my answer above.

I hope this clarifies further what is happening. Don't hesitate to send more questions.