Finding isolated points on a variety

[matthiashimmelmann]

Hey, I want to write a program that checks, whether a given point is (real) isolated. For example, given the variety V=Z(y(x-1), (x(x-1)) that is constituted by the union of (0,0) and the parametric line (1,y) inR^2, I look at the ring 0<x^2+y^2<=t for t>=0. By minimizing t, I suspect that I find the point on V\{(0,0)} closest to (0,0). If t converged to 0 during the optimization, that means that (0,0) is not isolated. Otherwise, for t>TOL, we conclude that (0,0) is indeed isolated, as no other point on the variety lies in close proximity. For this algorithm's implementation, I use the following code:

using SumOfSquares, DynamicPolynomials, MosekTools
@polyvar x[1:2]
S = @set x[1]*(x[2]-1) >= 0 && x[1]*(x[2]-1) <= 0 && x[2]*(x[2]-1) >= 0 && x[2]*(x[2]-1) <= 0 && x[2]^2+x[1]^2>=10^(-10)

model=SOSModel(Mosek.Optimizer)
@variable(model,t)
@constraint(model, con1, x[1]^2+x[2]^2<= t, domain = S)
@objective(model, Min, t )

optimize!(model)
@show termination_status(model)
@show value(t)

However, the model terminates with the code MathOptInterface.INFEASIBLE, so it does not actually determine, whether (0,0) is isolated. I would be very grateful if you helped me to figure out, how to get this model to find an optimal solution!

[blegat]

You want to find a t for which x[1]^2 + x[2]^2 >= t for all x such that x[1]^2 + x[2]^2 >= TOL and x belongs to the variety. Therefore, you need to write @constraint(model, x[1]^2 + x[2]^2 >= t, domain = S), you wrote >=. Also, the tolerance 1e-10 is maybe bit too low for a numerical SDP solver. I would suggest using 1e-1 to start with (since you know that it works for this example) and once it works you can try to see what is the limit for the precisions so that it does not run into numerical errors. Moreover, note that you can do @set x[1]*(x[2]-1) == 0 && x[2]*(x[2]-1) == 0 && x[2]^2+x[1]^2>=1e-1. So that we see that the basic semialgebraic set is defined partly by equalities and partly by inequalities and for the equality part, it will compute the Groebner basis check the equality between the polynomial and it's Sum-of-Squares decomposition in the reduced form. On other hand, for inequalities, we need to create a Sum-of-Squares polynomial for each inequality but we don't know which degree to use so it creates some concervatism. You'll probably need to play with the maxdegree keyword to increase the degree of the Sum-of-Squares polynomial used to for x[2]^2 + x[1]^2 - 1e-1 if you still get infeasibilities.

[matthiashimmelmann]

Thanks a lot for your help! The maxdegree keyword helped (curiously, 3 as input let the algorithm compute the optimal solution). Also, thanks for your remark concerning equality constraints involving Gröbner bases. Maybe I'll toy around with numerical approaches, if I look at very complex systems, when Gröbner bases become infeasible. Anyways, I'll try to push tolerance and complexity as far as I can.