Computing a point in a variety of a positive dimensional ideals

[DannyFeldman]

Problem Description

variety() supports only ideal of dimension zero, i.e., finite number of solutions.
It would be helpful if it would return at least one solution in such cases.
For example, in optimization problems, even argmin is an infinite set, we could conclude the minimum value.
It seems that commercial CAS such as Mathematica (FindInstance) or Maple (solve) support this feature.

Proposed Solution

The solution for existing software seems to be based on:
Patrizia Gianni, Properties of Gr¨obner bases under specialization, Lect. N. Comp. Sci. Berlin, Heidelberg, New York: Springer 378 (1987), 293–297. MR1033305 (91g:13032)

Michael Kalkbrenner, Solving systems of algebraic equations using Gr¨obner bases, Lect. N. Comp. Sci. Berlin, Heidelberg, New York: Springer 378 (1987), 282–292.

Alternatives Considered

An alternative solution might be to return a decomposition of the space of solutions. For example, "Algorithms in Real Algebraic Geometry" by Basu, Pollack and Roy suggests such a decomposition into constructible sets in Theorem 1.22.

Additional Information

No response

Is there an existing issue for this?

• I have searched the existing issues for a bug report that matches the one I want to file, without success.

[DannyFeldman]

@maxale, do you think variety() should return such a solution instead of "nothing"? Or should it be a different function such as "exampleSolution()"

[maxale]

@DannyFeldman: I think the default behavior of .variety() should not be changed, that is, it should not return nothing, but raise an error as it does now to notify the user about impossibility of computing the whole variety. As for the new functionality, it can be either kick in when an optional parameter to .variety() is provided (say, one_point=True), or come in a separate method such as .one_point(). I have no strong option here.

Also, if computing one point can be extended to two, three, etc., it may be worth it to provide a method like .points_iterator(n_points=+oo) that iterate over (some) points until n_points are generated.

[vneiger]

Agreeing with the above comment: this should not modify variety which has a clear description about returning the whole variety. My preference would go towards another method, e.g. an_element_of_variety (following the existing terminology an_element for many algebraic structures in Sage).

[vneiger]

Answering to #39475 (comment) , which I think belongs to this thread.

Thank you for the comments and information! in case of struggle with integration into SageMath, feel free to post the code somewhere and ask here on this issue if anyone would be able to help for the integration.

[DannyFeldman]

@vneiger Awesome! I will!

Another comment regarding motivation: our main application is minimizing (approximately) a polynomial function under polynomial constraints. Using binary search on the search space, the problem reduces to the decision problem "give me any solution to this system unless it has no solution (empty variety)". Many optimization problems in computer science, and computer vision can be reduced to these cases, where we do benchmark for minimizing a function, and less care about the "arg min" set which might be infinite.

[dimpase]

In the minimisation setting, there are results saying that you should look at the gradient ideal of the polynomial you minimise:
http://arxiv.org/pdf/math/0411342v2

In general, a lot depends on extra properties of your ideal $I$, or its Groebner basis (GB) you can get hold of.
E.g. if you have a lexicographic GB then you will have there a polynomial $p$ with $1+dim(I)$ variables, and then you can find a solution of $p=0$, and propagate it.

If your $I$ is homogeneous, intersecting its variety with a random linear subspace of codimension $dim(I)$ will give you a 0-dimensional ideal (add the linear forms forming a basis of this subspace to $I$).
This also says that you can homogenise $I$ and get a solution where some coordinates might be at infinity.

Also, in the approximate setting, there is a whole industry called "numerical algebraic geometry" which does just this - finding approximate solutions. See e.g. https://www3.nd.edu/~jhauenst/preprints/hsWhatIsNAG.pdf
Would be good if Sage got an interface to one of these methods, e.g. http://homepages.math.uic.edu/~jan/phcpy_doc_html/index.html (there was some work done in this direction, but I lost track of it somehow). @antonleykin

[DannyFeldman]

@dimpase ,

Wow, these are great research directions.
I investigated many polynomial solvers over the last year.
Including Resultant, Macaulay/Action matrix, and Homotopy.
There are gazilon of papers whose titles include the words "global optimum" "optimal" and "exact".
However, if you read each of them carefully, you find one of the two conditions:

1. There is no theorem for correctness, but only analysis that contains a hidden assumption very deep inside it. Say, that 0=1 "but our experimental results show that it is a reasonable assumption". Mainly in journals from applied communities such as computer vision.
2. There is a theorem for "global optimal solution", but only under some complicated conditions. It is usually not clear how to check these conditions, and this verificaiton problem is usually harder than the original problem.

there are results saying that you should look at the gradient ideal of the polynomial you minimise:
http://arxiv.org/pdf/math/0411342v2

Unfortunately, it seems that this reference (as most of the SOS papers) satisfies the second condition. Also, I am not sure that SOS supports complex polynomial system.

In general, a lot depends on extra properties of your ideal $I$, or its Groebner basis (GB) you can get hold of.

The problems I have in mind (mainly from computer vision) involves <5 variables, equations and degree. So, I really hope to find a general solution for such "simple" cases. The approximation factor can be calibrated or improved iteratively via the binary search.

If your $I$ is homogeneous...

Sounds really cool idea. Even when $I$ is non-homogenous, we can try to "homogonize" it, and then project. References (state-of-the-art or textbook for beginners) would be appreciated.

See e.g. https://www3.nd.edu/~jhauenst/preprints/hsWhatIsNAG.pdf

Sounds also very interesting. Although I am more suspiciuos. It sounds like Homotopy continuation (as in your last reference), which also satisfies the above conditions, as far as I know.

E.g. if you have a lexicographic GB then you will have there a polynomial $p$ with $1+dim(I)$ variables, and then you can find a solution of $p=0$, and propagate it.

This is exactly the solution we implemented, following the references in the previous comments. I will appreciate your comments on the pseudo code (coming soon).

[dimpase]

@DannyFeldman you need to say explicitly whether you are after any solution (that is, in an algebraic extension of the field) or only a real one. The latter are much harder - there the notion of the dimension of the ideal is different, and harder to work with.

Nothing I suggested is meant to address the real solutions specifically. For the real case, the reference would be Basu-Pollack-Roy book

The "hidden assumptions" are typical for the SOS-based methods, which have to work around the theoretical complexity issues with SDP solving (they want it to be polynomial-time, but nothing like this is known in general) to get a complexity bound.

On a practical side, msolve is the fastest of Groebner bases implementation available in SageMath (the interface is, as usual, not quite complete).

[antonleykin]

@DannyFeldman
In addition to all great pointers from @dimpase, there is a lot of recent work (including by yours truly) using some of these methods in computer vision. Pephaps e-mail me an example of a problem and a type of solution that you are thinking about --- I'd try to suggest where to look.

[dimpase]

E.g. if you have a lexicographic GB then you will have there a polynomial
....

This is exactly the solution we implemented, following the references in the previous comments. I will appreciate your comments on the pseudo code (coming soon).

Lexicographic is an overkill, though. An elimination order with appropriate block sizes will be much faster, in general.

[DannyFeldman]

@dimpase and @antonleykin, I am very excited by your comments and have many questions and answers. However, I will try to be focused for start and refer to more issues later.

Pephaps e-mail me an example of a problem
@antonleykin, many thanks. I will.

Lexicographic is an overkill,
I read this claim many times. However, all the solutions I know are based on the Elimination theory which is based on Lex Order. See the pseudo code in my following comment. Any ideas how to replace it would be appreciated.

you need to say explicitly whether you are after any solution (that is, in an algebraic extension of the field) or only a real one
Apologies, I thought that I mentioned that I am only interested in a single element in the variety of a complex polynomial system in C[x_1,..x_n]. The pseudo code and code are for these cases, but I guess it should work for any algebraically closed field.

For the real case, the reference would be [Basu-Pollack-Roy book]
In practice or SageMath, I guess this means computing CAD. However, CAD has reputation that is worse even than Groebner bases when it comes to running times. See next answer.

or only a real one. The latter are much harder
Exactly. While my real motivation is in the Euclidean space (machine learning and computer vision in particular), the plan is to assume that the dimension is zero, find all finite solutions and then return one of them. In fact, all the papers and books I see (SOS, homotopy, whatever) for solving polynomial systems assume finite number of solutions ("zero dimension" in the complex case")
If the dimension is more than one, I hope to reduce the system to zero dimension and hope that at least one real solution was left. For start, as a heuristic. Later, I hope to prove that it provably works, at least for specific cases.

[DannyFeldman]

poly_solve.pdf

Attached is a description of a pseudo code for computing a single element in a given basis for an idea in complex variables.
I intend to upload the SageMath after getting your comments.
First questions:

• Can we remove the Lex order by something more efficient, as @dimpase suggested?
• Is there a simple way to get "free parameters"? My post-doc Alex computed it "from scratch" in SageMath based on [1]. Maybe .dimension() of SageMath already compute them?
• Is there a nice way to test such algorithms? Currently we implemented kind of exhaustive search on random polynomial systems and evaluate the suggested solution.
• Lines 7--9 in Algorithm 1 just "change the order of variables". Is there a shorter but still formal way to write it?

[DannyFeldman]

The "hidden assumptions" are typical for the SOS-based methods, which have to work around the theoretical complexity issues >>with SDP solving (they want it to be polynomial-time, but nothing like this is known in general) to get a complexity bound.

Not at all. At least from my experience, and most of the papers I read. The running time of the SDP solvers is extremely fast in practice (and also in theory, after reasonable assumption on the input, such as conditioning). In fact, you can solve SOS without SDP:
https://arxiv.org/abs/1712.01792

The main problem I have, both in theory and practice, is that the optimality is guaranteed only when the hierarchy-level approaches infinity. Even in this case, it is only using many more assumptions. In practice, the action matrix explodes after two or three levels. Unlike other heuristics, the magics are that: 1) If and when it works -- you get a certificate of optimality, and 2) There are some very few but famous examples that solving an NP-hard problem in only the first levels gives a provable approximation in polynomial time. However, all the examples I know are not new and based on really hard tailoring of an existing approximation. See:
https://www.boazbarak.org/sos/prev/files/lec2d.pdf

But, at least for systems of a few (constant <10) polynomials, in a small degree and few variables, I expect Groebner basis to give all the solutions of a zero-dimensional system. I am still failing to do it with SageMath (will open a new "issue") but believe it will happen soon with your help.

[dimpase]

@DannyFeldman

Lexicographic is an overkill, I read this claim many times. However, all the solutions I know are based on the Elimination > theory which is based on Lex Order.

• Can we remove the Lex order by something more efficient, as @dimpase suggested?

yes, this is a standard procedure (not in Sage, though). Let $I\subset K[x_1,...x_k,y_1,...,y_m]:=T$ be an ideal and we want to compute $J=I\cap K[x_1,...x_k]:=S$. Any monomial order satisfying, for any $f\in T$ the property " $in(f)\in S$ implies $f\in S$", will give you a Groebner basis $G$ for $I$, so that $G\cap S$ is a Groebner basis for $J$. Cf. e.g. D. Eisenbud's "Commutative algebra; with a view towards algebraic geometry", Sect. 15.10.4.
Such orders are called elimination orders. These are arguably the true basis of Elimination Theory.
As implemented in Macaulay2: https://macaulay2.com/doc/Macaulay2/share/doc/Macaulay2/Macaulay2Doc/html/___Eliminate.html - SageMath has a (slow) interface to it, by the way.

[dimpase]

for systems of a few (constant <10) polynomials, in a small degree and few variables, I expect Groebner basis to give all the solutions of a zero-dimensional system. I am still failing to do it with SageMath (will open a new "issue") but believe it will happen soon with your help.

@DannyFeldman : Did you try using msolve backend? That's probably the best you can hope for in any current system.
Anyway, it's not clear what you mean by "solutions". It inevitably has to deal with using roots of a multivariate polynomial $f_0(t)$.
msolve can give you all solutions of a 0-dimensional system in form of $(\frac{f_1}{f_0^{\prime }},...,\frac{f_n}{f_0^{\prime }})$,
with $f_k$ univariate polynomials in $t$, which need to be evaluated at the roots of $f_0$. See Sect. 6 of https://msolve.lip6.fr/downloads/msolve-tutorial.pdf
Then you can plug numerically found roots of $f_0$ into these vectors of rational functions in $t$, to get all the roots of the system.
(IIRC there is no implementation in Sage to extract such solutions from msolve, but this is not hard to add)

[DannyFeldman]

@dimpase:

Did you try using msolve backend?

No, the student is still struggling with the compilation. Will keep you updated.

Anyway, it's not clear what you mean by "solutions"

Good point. In the context of this thread and the issue in the previous paragraph: any element in the variety in $C^n$ will do. Say, because we want to compute the global minimum of a given function, even if it has infinite such minimizers. Or using the binary search that I mentioned. The problem with my suggested pseudo code in the previous comment: the algorithm computes Groebner in the first lines of the code and sometimes hangs...

For example, , I would expect that groebner_basis()" would solve simple problem such as: Let $F,G\in R^{2\times 2}$. Minimize $||(F+y_{1}G)(x_1,x_2{)}^T+(0,y_2{)}^T|{|}^2$ over every pair of unit vectors $x,y\in R^2$.
I use Lagrange multipliers and then compute groebner_basis().

When $F$ and $G$ are given as "real" (say, random integers) it takes a second. However, when $F$ and $G$ are parameters, the computer hangs for days. In fact, the same happens in Mathematica or Maple. My suspicious is that the reason is that the system is easy to solve in the general case, when the dimension is zero (variety of finite size). However, in the de-generative cases, maybe when $F$ and $G$ have rank 0 or 1, there are infinite number of minimizers and the basis is huge. But the main motivation for implementing the pseudo code was to handle such cases. Maybe adding more constraints such as $t\ast det(F)-1=0$ would solve it.

Any suggestions would be appreciated.

[dimpase]

@dimpase:

Did you try using msolve backend?

No, the student is still struggling with the compilation. Will keep you updated.

Why is it a problem? I'd suggest to use Conda, that's the quickest way to get a development version of Sage ready. We'd be glad to help with this.

By the way, I just posted a PR to upgrade msolve to the latest version #39738

Anyway, it's not clear what you mean by "solutions"

Good point. In the context of this thread and the issue in the previous paragraph: any element in the variety in C n will do. Say, because we want to compute the global minimum of a given function, even if it has infinite such minimizers.

Do you have a function ${\mathbb{C}}^n\to \mathbb{R}$ to minimise? It's not very usual, that's
why I ask.

Or using the binary search that I mentioned. The problem with my suggested pseudo code in the previous comment: the algorithm computes Groebner in the first lines of the code and sometimes hangs...

For example, , I would expect that groebner_basis()" would solve simple problem such as: Let F , G ∈ R 2 × 2 . Minimize | | ( F + y 1 G ) ( x 1 , x 2 ) T + ( 0 , y 2 ) T | | 2 over every pair of unit vectors x , y ∈ R 2 . I use Lagrange multipliers and then compute groebner_basis().

I suppose you checked that it's well-posed, etc. Cause Lagrange multipliers need not always work.

When F and G are given as "real" (say, random integers) it takes a second. However, when F and G are parameters, the computer hangs for days. In fact, the same happens in Mathematica or Maple. My suspicious is that the reason is that the system is easy to solve in the general case, when the dimension is zero (variety of finite size). However, in the de-generative cases, maybe when F and G have rank 0 or 1, there are infinite number of minimizers and the basis is huge. But the main motivation for implementing the pseudo code was to handle such cases. Maybe adding more constraints such as t ∗ d e t ( F ) − 1 = 0 would solve it.

I'm not sure whether your example is a real one, or a toy, one, but just in case:
If $G$ is completely general then $F+y_{1}G$ may be replaced by unknown matrix $W$, so we'd minimise $|Wx+(0,y_2{)}^{\mathrm{\top }}{|}^2$ on the torus $|x{|}^2=|y{|}^2=1$, and it's obvious that at a global minimum you must have $y_2=0$. This leaves you with the classical problem of minimizing $|Wx|$ on the unit circle, i.e. finding the minimal eigenvalue of $W$, and you can write down the answer by hand.
PS. the hard case is when $W$ has repeated eigenvalues - something one can fix by stipulating non-vanishing of the discriminant of the char. poly. of $W$.

[dimpase]

The "hidden assumptions" are typical for the SOS-based methods, which have to work around the theoretical complexity issues >>with SDP solving (they want it to be polynomial-time, but nothing like this is known in general) to get a complexity bound.

Not at all. At least from my experience, and most of the papers I read. The running time of the SDP solvers is extremely fast in practice (and also in theory, after reasonable assumption on the input, such as conditioning). In fact, you can solve SOS without SDP: https://arxiv.org/abs/1712.01792

In my research experience even solving LPs in "practice", e.g. while computing Delsarte bounds on binary code sizes, can be impossible (without a special arbitrary precision solvers)
with relatively few, ~100, variables. Indeed, conditioning. Similar kind of SDP-based bounds are even worse. You need arbitrary precision, but naively it blows up, memory-wise, after few iterations.
Things like certifying that a positive polynomial is SOS also often blow up as soon as your polynomial is close to the boundary of the SOS cone in the bigger cone of nonnegative polynomials.

[DannyFeldman]

In my research experience even solving LPs in "practice"
Right, it might be related to the "initialization seed" issues. There are probably other many reasons to to use SOS for such simple and small problems.

This is why I am so surprised that Groebner() gets lost already there. Hope that at least msolve would change it. Meanwhile I try to implement a solution or at least to understand why is it so hard. It must be related to the degenerate cases, but I am not sure why they cannot be removed in advance.

[DannyFeldman]

@dimpase , thank you for the detailed reply.

Why is it a problem?

Not sure. Alex writes a nice report for installation a complete version of Sagemath where you can actually debug the groebner_basis function (there were some bugs in the latest version). Will update with a summary that hopefully help others.

Do you have a function C n → R to minimise? It's not very usual, that's why I ask.
I actually have a function from R^4 to R as below. The idea is to use complex (which are simpler due to Groebner, no CAD). This is by using Lagrange multipliers, and hope to get a variety of a finite size (and filter the real ones), or just use the binary search that i mentioned above and hope that the algorithm in the pseudo code that I posted would return a real and not complex item in the variery.

I suppose you checked that it's well-posed, etc. Cause Lagrange multipliers need not always work.

At least in the general (and interesting case), there are finite number of extremum and Lagrange seems ok.

I'm not sure whether your example is a real one,

It is a real problem but I am actually looking for a generic solution for future problems, not a trick. Here is the motivation:
The input to the "noisy Perspective-3-Points" is a set of 3 points in R^3, paired with 3 lines that intersect the origin. The goal is keep the first points on the first line, the second point on the second line, and minimize the distance from the third point to the third line. It is (not so) easy to verify that the set of all possible distances is
{\norm{(F+y_1G)x +(0,cy_2)^T}\mid x,y\in S}, where $S$ is the unit sphere in $R^2$, $c\ge 0$ and $F,G$ are matrices that depend on the input.

Variants to the simpler 2d case y=(0,1)^T can be found here:
https://arxiv.org/abs/1807.08446

If G is completely general then F + y 1 G may be replaced by unknown matrix W ,
Firstly, F and G are the input parameters (of a given instance) and not part of the optimization problem (on x and y). So I can't see how you can assume this. I probably miss something.
Secondly, at least in the real case that I have in mind, y_1=<1 due to the constraint that y is a unit vector.

This leaves you with the classical problem of minimizing | W x | on the unit circle

This is too good to be true, so if you are right (see previous comment) I probably missed something in the problem. Since even when the smallest norm (distance) above is zero, the problem is the (noise-free) P3P which reduces to solving a single univariate quartic polynomial.

[DannyFeldman]

@vneiger,

in case of struggle with integration into SageMath, feel free to post the code somewhere and ask here on this issue if anyone would be able to help for the integration.

Thank you so much for the willing to help. Dr. Verner is currently working on installing (the full version of) SageMath, and would appreciate your help. See:
[(https://github.com//discussions/39272#discussioncomment-12556848)[

The goals are to: 1) add debug (or timer) info in groebner_basis, 2) Run msolve, 3) Implement the pseudo code that I uploaded above.

[dimpase]

please refrain from discussing at the discussion there, it on a different, concrete, topic! We teach our 1st year undergrads to not hijack existing topics like this :-)

I've opened a new issue, #39740, to deal with this.

[DannyFeldman]

please refrain from discussing at the discussion there, it on a different, concrete, topic! We teach our 1st year undergrads to not hijack existing topics like this :-)

So sorry about this. Alex opened it but we are both new in this business. Will be more careful next time.

I've opened a new issue, #39740, to deal with this.

Many thanks. Fingers crossed...

[DannyFeldman]

Dear Dima, I apologize for bothering you with this again and my ignorance. Can you please confirm that the place for the question of Alex is: #39740 Alex, I guess that Dima meant that you posted the question in the wrong (existing) case. Not the wrong forum. Your report and tutorial for the installation might also be used by the community (after verification). Thank you both, Danny

…

On Thu, 20 Mar 2025 at 11:41, Alexander Verner ***@***.***> wrote: Yes, it appears here: #39740. But this contradicts his disapproval of posting in the GitHub forum. чт, 20 мар. 2025 г. в 11:37, Dan Feldman ***@***.***>: > Hi, Alex > Dima wrote below that he "opened a new issue to deal with this." > Maybe this is why your new cases are being deleted. > > Do you find this new open case? > > ---------- Forwarded message --------- > From: Dima Pasechnik ***@***.***> > Date: Thu, 20 Mar 2025 at 00:09 > Subject: Re: [sagemath/sage] Computing a point in a variety of a positive > dimensional ideals (Issue #39682) > To: sagemath/sage ***@***.***> > Cc: DannyFeldman ***@***.***>, Mention < > ***@***.***> > > > please refrain from discussing at the discussion there, it on a > different, concrete, topic! We teach our 1st year undergrads to not hijack > existing topics like this :-) > > I've opened a new issue to deal with this. > > — > Reply to this email directly, view it on GitHub > <#39682 (comment)>, > or unsubscribe > <https://github.com/notifications/unsubscribe-auth/ALAQHRMGD3ZAIDPGZK4JPVT2VHTHTAVCNFSM6AAAAABY35YZZWVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDOMZYGI4DGOBZG4> > . > You are receiving this because you were mentioned.Message ID: > ***@***.***> > [image: dimpase]*dimpase* left a comment (sagemath/sage#39682) > <#39682 (comment)> > > please refrain from discussing at the discussion there, it on a > different, concrete, topic! We teach our 1st year undergrads to not hijack > existing topics like this :-) > > I've opened a new issue to deal with this. > > — > Reply to this email directly, view it on GitHub > <#39682 (comment)>, > or unsubscribe > <https://github.com/notifications/unsubscribe-auth/ALAQHRMGD3ZAIDPGZK4JPVT2VHTHTAVCNFSM6AAAAABY35YZZWVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDOMZYGI4DGOBZG4> > . > You are receiving this because you were mentioned.Message ID: > ***@***.***> >

[dimpase]

@DannyFeldman no worries, I actually opened #39740 myself - copying Alex's question there. Perhaps the title should be changed the, to reflect that there is also a building issue.

[vneiger]

@vneiger,

in case of struggle with integration into SageMath, feel free to post the code somewhere and ask here on this issue if anyone would be able to help for the integration.

Thank you so much for the willing to help. Dr. Verner is currently working on installing (the full version of) SageMath, and would appreciate your help. See: [(https://github.com//discussions/39272#discussioncomment-12556848)[

The goals are to: 1) add debug (or timer) info in groebner_basis, 2) Run msolve, 3) Implement the pseudo code that I uploaded above.

Sorry for not answering earlier, it was very busy days. I see some help has been provided already, thank you @dimpase ; I'm subscribing to this other issue to keep informed if help is needed.