sagemath / sage

Main repository of SageMath
https://www.sagemath.org
Other
1.44k stars 480 forks source link

multivariate factorization over QQbar #25390

Closed BrentBaccala closed 5 years ago

BrentBaccala commented 6 years ago

We can currently (#8544) factor univariate polynomials over QQbar:

sage: R.<x>=QQbar[]
sage: (x^2+1).factor()
(x - I) * (x + I)

We can factor over multivariate rings, if the polynomial is actually univariate:

sage: R.<x,y>=QQbar[]
sage: (x^2+1).factor()
(x - I) * (x + I)

But we can't do full "absolute factorization" (that's what it's called in the literature):

sage: R.<x,y>=QQbar[]
sage: (x^2+y^2).factor()
...
NotImplementedError: proof = True factorization not implemented.  Call factor with proof=False.
sage: (x^2+y^2).factor(proof=False)
...
TypeError: no conversion of this ring to a Singular ring defined

This ticket implements absolute factorization by calling Singular's absfact.lib, which works for rings over Q (or transcendental extensions thereof), and uses a technique described in the attached document to handle number fields.

I flagged it major because polynomial factorization is a fundamental feature that blocks the implementation of other things, like primary decomposition of ideals.

Depends on #25351

CC: @tom111

Component: algebra

Author: Brent Baccala

Branch/Commit: 91d467e

Reviewer: Travis Scrimshaw

Issue created by migration from https://trac.sagemath.org/ticket/25390

BrentBaccala commented 6 years ago

Commit: 27f109d

BrentBaccala commented 6 years ago

Branch: public/25390

BrentBaccala commented 6 years ago
comment:1

I added code to call Singular's absFactorize if the polynomial's coefficients are in QQ.

Here's what's now possible:

sage: R.<x,y> = QQbar[]

sage: L = (x^2+y^2).factor()
sage: L
((-1*I)*x - y) * (1*I*x - y)
sage: L.value()
x^2 + y^2

sage: p = (-7*x^2 + 2*x*y^2 + 6*x + y^4 + 14*y^2 + 47)*(5*x^2+y^2)^3*(x-y)^4
sage: F = p.factor()
sage: F
...
sage: F.value() == p
True

But something like this still doesn't work (because the coefficients aren't in QQ):

sage: (x+QQbar(sqrt(2))*y).factor()

It's something, which is better than nothing, so I think we should move it towards master, but there's still a lot of work to be done before this ticket is closed.


New commits:

27f109dTrac #25390: use Singular 'absfact' to factor multivariate polynomials over QQ
BrentBaccala commented 6 years ago

Author: Brent Baccala

tscrim commented 6 years ago
comment:2

I agree that small improvements are good. We can always do followup tickets to add more functionality. I do have some comments:

Can we use libsingular instead of calls to singular (which uses the pexpect interface and is much slower)?

Doc tweaks:

         ALGORITHM:

-        Uses Singular's `absfact` library.
+        Uses Singular's ``absfact`` library.

-        TODO:
+        .. TODO::

-        Implement absolute factorization over number fields
+            Implement absolute factorization over number fields.

Also if you could split the long output line so that it is at most ~80 chars per line.

Instead of elem_map, you could just use elem_dict.__getitem__.

I don't see the need for the internal functions polynomial_map and reverse_polynomial_map. These little trivial functions make the code harder to follow.

In Python, usually asserts are with spaces, e.g., assert minpoly.degree() == 0.

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

e6124e2Trac #25390: convoluted code that I want to save before cleaning
b34f2b3Trac #25390: correctly handle factorization over number fields
2879cc0Trac #25390: bug fix (and add a test case)
747861eTrac #25390: code cleanup suggested by tscrim's review
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Changed commit from 27f109d to 747861e

BrentBaccala commented 6 years ago
comment:5

Over the weekend I remembered a technique that I read years ago for factoring over a number field - factor the norm of the polynomial, which has rational coefficients and is a multiple of the original polynomial.

So, this ticket is moving towards closed more quickly that I had thought possible. I think it can now factor any polynomial over QQbar.

The algorithm is tricky enough that I'm planning to write up a short paper describing it, post it on arxiv, and link to it from the documentation. If anybody has a better suggestion for how to incorporate a three page paper explaining the algorithm of a Sage method, please let me know.

I took most of tscrim's suggestions, but I can't figure out how to access a Singular variable from libsingular. I posted a question on asksage, but so far no answer.

I also intend to implement factorization over AA before flipping this ticket to needs_review.

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Changed commit from 747861e to 247f944

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

b2657e2Trac #25390: simplify factorization logic
580c11eTrac #25390: enhance multivariate factorization to handle AA as well as QQbar
be078acTrac #25390: ensure that multivariate factorization factors are monic
7b6c93bTrac #25390: add a test case
247f944Trac #25390: slight code rearrangement
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

efdaac1Merge tag '8.3.beta5' into public/25390
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Changed commit from 247f944 to efdaac1

BrentBaccala commented 6 years ago
comment:8

Still haven't figured out how to access the Singular variable from libsingular; got an answer to that SageMath question that made sense, but nbruin (the responder) said that it was "tricky" and he was right.

I spent two weeks trying to figure out to make this code work without doing division in polynomial rings over QQbar but it's a lot harder than I thought. I'm attaching a PDF that describes the progress I've made and the problems that I've had.

Right now, the code is fairly straightforward, but does require polynomial division over QQbar, and therefore requires Trac #25351 to be applied before it works. I've spent enough time trying to break that dependency, so I'm flipping this ticket to needs_review and leaving it the way it is. See the attached PDF for more details.

BrentBaccala commented 6 years ago

Dependencies: #25351

tscrim commented 6 years ago
comment:9

Attachment: qqbar.pdf.gz

I didn't realize that my suggestion would have been so complicated. Thank you for looking into it. I have two small additional nitpicks, but otherwise this ticket is a positive review assuming #25351 (which I cannot review as that is much more involved mathematically, you should cc some people who work more in that area, such as jdemeyer).

You do not need to create a list:

norm_f = prod([numfield_f.map_coefficients(h) for h in numfield.embeddings(QQbar)]).change_ring(QQ)
norm_f = prod(numfield_f.map_coefficients(h) for h in numfield.embeddings(QQbar)).change_ring(QQ)
-        REFERENCE::
+        REFERENCES:

-            Geddes, et. al, "Algorithms for Computer Algebra", Section 8.8
+        - Geddes, et. al, "Algorithms for Computer Algebra", Section 8.8

although I think it would be better to reference the book in the master ref file and use - [Geddes1234]_ Section 8.8 as the reference.

nbruin commented 6 years ago
comment:10

If at all possible, you'll probably want to avoid computing norms by multiplying complex embeddings together. Especially when you have your polynomial over a number field, it's probably much better to avoid floats altogether. For instance, if your number field is given as K=QQ[t]/(h(t)) (with h(t) monic) and you have a polynomial f(x) in K[x], then you can represent your polynomial as F(x,t) in QQ[x,t], via the isomorphism K[x] ~ Q[x,t]/(h(t)). Then you have that

Norm_(K[x]/Q[x]) (f) = Resultant(F,h,t)

Code for computing resultants is usually quite optimized. You can experiment a bit to see which approach works best, but I'd expect the resultant to work quite well.

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Changed commit from efdaac1 to cbc51c4

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

3adbd36Trac #25390: use resultant to compute polynomial norm
9d6353aTrac #25390: bug fix (if polynomial names conflict with number field names)
cbc51c4Trac #25390: update references
BrentBaccala commented 6 years ago
comment:12

I made the changes suggested by tscrim and nbruin; thanks for your feedback.

Using the resultant was more of a pain that I thought it should be. The difficulty lay mainly in needing to introduce a new variable, different from the others. If anybody can suggest a better way that what I came up with, please let me know.

nbruin commented 6 years ago
comment:13

Variable names might be haunting us here:

sage: k=NumberField(x^2+1,"x")
sage: R=PolynomialRing(k,"x,y")
sage: P=PolynomialRing(QQ,k.gens()+R.gens())
ValueError: variable name 'x' appears more than once

Unless you can guarantee the variable names of the number field and the polynomial ring aren't clashing, it seems you cannot go about the problem this way. While it would be reasonable to assume that a user wouldn't intentionally throw this kind of awful naming at you, it's quite conceivable that automatic code would end up doing something like this. Example:

sage: matrix(k,2,2,[k.0,0,0,1]).charpoly()
x^2 + (-x - 1)*x + x

So it would seem the names of your variables for the trivariate ring should be unrelated to the names that occur for the others. Assuming you have a polynomial g in a ring P over a number field k over QQ, I think you can do your conversion with something along the lines of

R=g.parent()
k=R.base_ring()
P=PolynomialRing(QQ,len(k.gens()+R.gens()),'T')   
G=sum(P({(0,)+tuple(k):1})*v.polynomial()(P.0) for k,v in g.dict().iteritems())
NG=G.resultant(k.polynomial()(P.0),P.0)
RQ=PolynomialRing(QQ,R.gens())
NGQ=NG((0,)+RQ.gens())

Doing your conversions like this has the advantage that you don't rely on names of generators at all.

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

5b23d3fTrac #25390: improve norm calculation code
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Changed commit from cbc51c4 to 5b23d3f

BrentBaccala commented 6 years ago
comment:15

OK, I re-coded it the way nbruin suggested.

nbruin commented 6 years ago
comment:16
+            norm_f = sum(norm_ring({tuple(k)[1:]:v})
+                         for k,v in norm_flat.dict().iteritems())

Please don't think that unpacking a polynomial into a dictionary is the best way of handling it! To "lift" number field elements, I didn't see another way, but in this case

norm_f = norm_flat((0,)+norm_ring.gens())

would work just fine, as far as I can see. If you have a good reason to use the dict upacking (is it faster?) then it's fine, of course; but otherwise using the simpler code is probably preferable.

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

a3eb8efTrac #25390: simplify norm calculation code
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 6 years ago

Changed commit from 5b23d3f to a3eb8ef

BrentBaccala commented 6 years ago
comment:18

I thought there had to a better way to do it!

mkoeppe commented 6 years ago
comment:19

Attachment: qqbar.tex.gz

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 5 years ago

Changed commit from a3eb8ef to ac30e3a

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 5 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

ac30e3aMerge tag '8.4' into public/25390
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 5 years ago

Changed commit from ac30e3a to 94a6c6c

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 5 years ago

Branch pushed to git repo; I updated commit sha1. Last 10 new commits:

730a40dTrac #25351: remove @handle_AA_and_QQbar decorator from reduce method
a9fde71Trac #25351: remove @handle_AA_and_QQbar decorator from _normal_basis_libsingular
4f547afTrac #25351: Put @handle_AA_and_QQbar decorator on groebner_basis(), and remove
d187260Trac #25351: improve @handle_AA_and_QQbar to handle keyword arguments
f83f2bdTrac #25351: add test cases to methods decorated with @handle_AA_and_QQbar
fb145a9Trac #25351: modify a test case to check keywords arguments to @handle_AA_and_QQbar
758e4eaMerge tag '8.7.rc0' into public/25351
7ba975eMerge tag '8.7' into public/25351
bd8ab1bMerge tag '8.8.beta2' into public/25351
94a6c6cMerge branch 'public/25351' into public/25390
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 5 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

1a3b28cMerge tag '8.8.beta7' into public/25390
7830854Trac #25351: ensure Python 3 compatibility
ebfbf2dMerge tag '8.8.beta5' into public/25351
b602cb1Merge tag '8.8.beta7' into public/25351
ba45dd6Merge branch 'public/25351' into public/25390
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 5 years ago

Changed commit from 94a6c6c to ba45dd6

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 5 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

10554b0Trac #25390: fix some doctests
f36b116Merge tag '8.8.rc0' into public/25390
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 5 years ago

Changed commit from ba45dd6 to f36b116

BrentBaccala commented 5 years ago

Description changed:

--- 
+++ 
@@ -26,8 +26,6 @@
 TypeError: no conversion of this ring to a Singular ring defined

-I was going to implement this using Singular's absfact.lib, but then I realized that those routines only work for rings over Q (or transcendental extensions thereof), and not over number fields.

-So, it's not just as easy as calling absFactorize. +This ticket implements absolute factorization by calling Singular's absfact.lib, which works for rings over Q (or transcendental extensions thereof), and uses a technique described in the attached document to handle number fields.

I flagged it major because polynomial factorization is a fundamental feature that blocks the implementation of other things, like primary decomposition of ideals.

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 5 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

64dd01eMerge tag '8.9.rc1' into public/25390
4a955fbTrac #25390: fix failing doctest
3fd38a8Trac #25390: remove a duplicate check for _factor_multivariate_polynomial
dd4840cTrac #25390: add a TODO item - maybe use univariate factorization code
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 5 years ago

Changed commit from f36b116 to dd4840c

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 5 years ago

Changed commit from dd4840c to a9c1e8b

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 5 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

a9c1e8bTrac #25390: replace iteritems() with items() for Python 3 compatibility
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 5 years ago

Changed commit from a9c1e8b to 6305cc1

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 5 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

6305cc1Merge tag '9.0.beta3' into public/25390
tscrim commented 5 years ago
comment:28

Two minor things: The first is the .. TODO:: block should be indented. The second is more stylistic, but I would avoid the blankline after the start of the for loops.

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 5 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

91d467eTrac #25390: indent TODO block and remove some blank lines
7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 5 years ago

Changed commit from 6305cc1 to 91d467e

BrentBaccala commented 5 years ago
comment:30

Replying to @tscrim:

Two minor things: The first is the .. TODO:: block should be indented. The second is more stylistic, but I would avoid the blankline after the start of the for loops.

OK, done.

I didn't remove the blank lines on the second for loop because it's so long, I feel the blank lines improve it readability.

tscrim commented 5 years ago

Reviewer: Travis Scrimshaw

tscrim commented 5 years ago
comment:31

Thanks. I disagree about the space because I have never really seen that elsewhere in Python code and it looks strange to be because of that. Yet, I don't think it is important enough to hold up the ticket.

vbraun commented 5 years ago

Changed branch from public/25390 to 91d467e