search

oscar-system / Oscar.jl

A comprehensive open source computer algebra system for computations in algebra, geometry, and number theory.
https://www.oscar-system.org
Other
338 stars 120 forks source link

Difficulties with radical computation #3711

Open HechtiDerLachs opened 4 months ago

HechtiDerLachs commented 4 months ago

Do I = load("non_radical_ideal.txt") with this file.

For me radical(I) does not finish, even though, in my opinion, the ideal is rather tame.

Interestingly, the following works:


julia> dec = primary_decomposition(I);

julia> intersect([p for (p, _) in dec]...)
Ideal generated by
  (y//z)
  (x//z)
  θ_2^2*θ_1^2 - θ_2^2 + θ_1^2 + 1
  2*θ_2^4 - 4*θ_2^2 + θ_1^4 + 2*θ_1^2 + 3
  2*θ_2^4*t - 4*θ_2^2*t + θ_1^4*t + 2*θ_1^2*t + 3*t
  θ_2^2*θ_1^2*t^2 - θ_2^2*t^2 + θ_1^2*t^2 + t^2
  -4*θ_2^2 + θ_1^6 + θ_1^4 + 7*θ_1^2 + 3
  θ_2^2*θ_1^2*t^3 - θ_2^2*t^3 + θ_1^2*t^3 + t^3
  2*θ_2^4*t^3 - 4*θ_2^2*t^3 + θ_1^4*t^3 + 2*θ_1^2*t^3 + 3*t^3
  t^8 + 1
  -4*θ_2^2*t^2 + θ_1^6*t^2 + θ_1^4*t^2 + 7*θ_1^2*t^2 + 3*t^2
  -4*θ_2^3*t + θ_2*θ_1^6*t + θ_2*θ_1^4*t + 7*θ_2*θ_1^2*t + 3*θ_2*t
  2*θ_2^5*t^4 - 4*θ_2^3*t^4 + θ_2*θ_1^4*t^4 + 2*θ_2*θ_1^2*t^4 + 3*θ_2*t^4
  -4*θ_2^2*t^3 + θ_1^6*t^3 + θ_1^4*t^3 + 7*θ_1^2*t^3 + 3*t^3
  2*θ_2^6*t^3 - 2*θ_2^2*t^3 + θ_1^4*t^3 + 3*t^3
  2*θ_2^4*t^6 - 4*θ_2^2*t^6 + θ_1^4*t^6 + 2*θ_1^2*t^6 + 3*t^6
  -4*θ_2^2*t^4 + θ_1^6*t^4 + θ_1^4*t^4 + 7*θ_1^2*t^4 + 3*t^4
  -4*θ_2^3*t^3 + θ_2*θ_1^6*t^3 + θ_2*θ_1^4*t^3 + 7*θ_2*θ_1^2*t^3 + 3*θ_2*t^3

So computing the radical as the intersection of associated primes works. This seems wrong, doesn't it? Also @simonbrandhorst has reported difficulties when computing radicals today. Has there been some regression?

@hannes14 ? @wdecker ?

HechtiDerLachs commented 4 months ago

I'm wondering whether there is some deeper issue with calls to Singular routines. Also the timeout of the tests in #3705 suggests this. Plus I'm sure that the latter went through on my machine before I updated to the latest master today.

wdecker commented 4 months ago

I cannot reproduce this result:

julia> @time radical(I);
  0.000009 seconds

julia> @time is_radical(I)
  0.000012 seconds
true
fieker commented 4 months ago

On Wed, May 15, 2024 at 06:35:22AM -0700, Wolfram Decker wrote:

I cannot reproduce this result:

julia> @time radical(I); 0.000009 seconds

julia> @time is_radical(I) 0.000012 seconds true

The radical is cached on the ideal, hence immedaite on 2nd call It takes 12sec for me to get it

-- Reply to this email directly or view it on GitHub: https://github.com/oscar-system/Oscar.jl/issues/3711#issuecomment-2112552867 You are receiving this because you are subscribed to this thread.

Message ID: @.***>

HechtiDerLachs commented 4 months ago

I just tried this again on latest master and still the computation does not finish:

R, (x,y,z,v,w) = QQ[:x, :y, :z, :v, :w]
g = QQMPolyRingElem[z^3-z*w^8-z - v^2, 3*z^2-w^8-1, v, z*w, x^8+1, -x^6+x^4-x^2-y^2]
I = ideal(R, g)
@time radical(I)

I'm running on Ubuntu, julia version 1.10.0, commit 5aa97ef175df57807ad1906e5b4410ea6f74d6e4 on master.

HechtiDerLachs commented 4 months ago

Update: I investigated this a bit with @fieker and it seems that on my machine julia decides to delegate this to one of many further processes. Eventually the computation finishes with 

julia> @time radical(I)
205.120327 seconds (1.44 k allocations: 36.656 KiB)
Ideal generated by
  v
  z
  x^2*y^2 + x^2 + y^2 - 1
  2*x^4 - 4*x^2 + y^4 - 2*y^2 + 3
  4*x^2 + y^6 - y^4 + 7*y^2 - 3
  w^8 + 1

It seems that almost no memory is used for this computation despite the long duration. This might be due to allocations in other processes not being counted?

Further version info on my machine:

(@v1.10) pkg> st
Status `~/.julia/environments/v1.10/Project.toml`
  [e30172f5] Documenter v1.4.1
  [3e1990a7] Hecke v0.31.5
  [f1435218] Oscar v1.1.0-DEV 
  [c46f51b8] ProfileView v1.7.2
  [295af30f] Revise v3.5.14

In Singular the same computation takes only fractions of a second.

joschmitt commented 4 months ago

I can confirm these timings. (The OSCAR ones.)

thofma commented 4 months ago

I cannot reproduce on macOS. Just as a breadcrumb, when it hangs on Linux, it seems to be hanging in

signal (10): User defined signal 1
__select at /lib/x86_64-linux-gnu/libc.so.6 (unknown line)
_Z12slStatusSsiLP6slistsi at /home/basdsad/.julia/artifacts/d742f40bb63b9380341851126ab31d8d64d59941/lib/libSingular.
so (unknown line)

This would be

slStatusSsiL(slists*, int)

but I have no clue about Singular internals.

thofma commented 4 months ago

Next breadcrumb: I cannot reproduce the error with Oscar 1.0.2. There it behaves normal. We could produce a Singular.jl only example and then do a bisect.

HechtiDerLachs commented 4 months ago

We have discussed this over here in KL a bit. I did not understand the issue 100%, but @hannes14 said something that certain library routines try to fork extra processes and do computations in parallel. Doing this from inside Oscar actually spawns various new julia processes and, as it seems, the communication to these processes is not handled correctly. So it is a game of luck whether you get a result at all and, according to @fieker, if you do, then most of the time is spent on waiting for other processes to finish.

fingolfin commented 4 months ago

@joschmitt what do you mean with "The OSCAR ones" ?

On macOS with current OSCAR master 3927f4bae77df98791af3517905b4b3bcd115bbf, like for @thofma, all of this is fast (on an M1 MacBook Pro):

julia> R, (x,y,z,v,w) = QQ[:x, :y, :z, :v, :w];

julia> g = QQMPolyRingElem[z^3-z*w^8-z - v^2, 3*z^2-w^8-1, v, z*w, x^8+1, -x^6+x^4-x^2-y^2];

julia> I = ideal(R, g);

julia> @time radical(I);
  3.340391 seconds (1.37 k allocations: 34.805 KiB)

Running this on tutulla (a Linux server with AMD CPU in Kaiserslautern), it indeed seems to take minutes...

But here is a funny thing: when I abort this by pressing Ctrl-C a couple times, and then (in the same session) enter all inputs again, it is fast, and stays fast! So I now got this:

julia> R, (x,y,z,v,w) = QQ[:x, :y, :z, :v, :w];

julia> g = QQMPolyRingElem[z^3-z*w^8-z - v^2, 3*z^2-w^8-1, v, z*w, x^8+1, -x^6+x^4-x^2-y^2];

julia> I = ideal(R, g);

julia> @time radical(I);
  3.482706 seconds (264.49 k allocations: 18.174 MiB, 5.77% compilation time)
fingolfin commented 4 months ago

@hannes14 @HechtiDerLachs so can we tell Singular to not try to spawn additional processes?

HechtiDerLachs commented 4 months ago

@hannes14 agreed to look into this. But as far as I understand, there is no simple global flag for that and you really need to adjust the actual library code for that.

joschmitt commented 4 months ago

@joschmitt what do you mean with "The OSCAR ones" ?

I meant that it also takes ~200 seconds for me with OSCAR master (and I did not try Singular).

thofma commented 4 months ago

What about the previous Singular/Singular.jl version, where it just works? Is it just luck, that the example works with Oscar 1.0.2 or did something change?

Edit: Nevermind, it has the same problems I think.

HechtiDerLachs commented 4 months ago

I am wondering whether the following call to is_prime might be affected by the same issue:

Load this polynomial into Oscar via

f = load("f.txt")
I = ideal(parent(f), f)
is_prime(I)

The computation does not finish for me. My current versions are

(@v1.10) pkg> st
Status `~/.julia/environments/v1.10/Project.toml`
  [c3fe647b] AbstractAlgebra v0.41.6
  [e30172f5] Documenter v1.4.1
  [3e1990a7] Hecke v0.32.0 `../../../Hecke.jl`
  [f1435218] Oscar v1.1.0-DEV `../../../Oscar.jl`
  [c46f51b8] ProfileView v1.7.2
  [295af30f] Revise v3.5.14

The problem did not occur before my latest updates of Hecke and AA today, i.e. the same code was running without difficulties before. I'm on a custom branched which is rebased on today's master branch.

JohnAAbbott commented 4 months ago 
R, (x,y,z,v,w) = QQ[:x, :y, :z, :v, :w]
g = QQMPolyRingElem[z^3-z*w^8-z - v^2, 3*z^2-w^8-1, v, z*w, x^8+1, -x^6+x^4-x^2-y^2]
I = ideal(R, g)
@time radical(I)

The same computation using CoCoA took less than 1s... unless I have misunderstood.

hannes14 commented 4 months ago

Singular also can compute it in less than 1 s, Also tried in a loop (so it does not depend on a "bad random" choice).

fieker commented 4 months ago

On Wed, May 22, 2024 at 04:59:36AM -0700, JohnAAbbott wrote:

R, (x,y,z,v,w) = QQ[:x, :y, :z, :v, :w]
g = QQMPolyRingElem[z^3-z*w^8-z - v^2, 3*z^2-w^8-1, v, z*w, x^8+1, -x^6+x^4-x^2-y^2]
I = ideal(R, g)
@time radical(I)

The same computation using CoCoA took less than 1s... unless I have misunderstood.

It's trying to do stuff in parallel again -- Reply to this email directly or view it on GitHub: https://github.com/oscar-system/Oscar.jl/issues/3711#issuecomment-2124617048 You are receiving this because you were mentioned.

Message ID: @.***>

fingolfin commented 4 months ago

@hannes14 I just tried this with the new Singular.jl. On my macBook, I now got a crash after trying a bunch of radical computations (though to be fair, this issue may have existed before, I just didn't try it ):

julia> R, (x,y,z,v,w) = QQ[:x, :y, :z, :v, :w]
(Multivariate polynomial ring in 5 variables over QQ, QQMPolyRingElem[x, y, z, v, w])

julia> g = QQMPolyRingElem[z^3-z*w^8-z - v^2, 3*z^2-w^8-1, v, z*w, x^8+1, -x^6+x^4-x^2-y^2]

6-element Vector{QQMPolyRingElem}:
 z^3 - z*w^8 - z - v^2
 3*z^2 - w^8 - 1
 v
 z*w
 x^8 + 1
 -x^6 + x^4 - x^2 - y^2

julia> @time radical(ideal(R, g))
  3.644767 seconds (1.39 k allocations: 36.055 KiB)
Ideal generated by
  v
  z
  x^2*y^2 + x^2 + y^2 - 1
  2*x^4 - 4*x^2 + y^4 - 2*y^2 + 3
  4*x^2 + y^6 - y^4 + 7*y^2 - 3
  w^8 + 1

julia> @time radical(ideal(R, g))
  3.412770 seconds (1.39 k allocations: 36.055 KiB)
Ideal generated by
  v
  z
  x^2*y^2 + x^2 + y^2 - 1
  2*x^4 - 4*x^2 + y^4 - 2*y^2 + 3
  4*x^2 + y^6 - y^4 + 7*y^2 - 3
  w^8 + 1

julia> @time radical(ideal(R, g))
// ** killing the basering for level 0

[88186] signal (11.2): Segmentation fault: 11
in expression starting at REPL[33]:1
_Z15ii_CallLibProcMPKcPPvPiP8ip_sringRi at /Users/mhorn/.julia/artifacts/26aefedef93010c5737ba6d78fac31dfa46979b9/lib/libSingular-4.4.0.dylib (unknown line)
Allocations: 39013484 (Pool: 38952023; Big: 61461); GC: 562
/Users/mhorn/.local/bin/j10: line 2: 88186 Segmentation fault: 11  julia +1.10 --proj=p/1.10 "$@"

On Linux, I tried on tutulla, and... the original issue doesn't seem to be fixed? Or at least not fully, as it is slower by a factor 10 compared to my MacBook

julia> R, (x,y,z,v,w) = QQ[:x, :y, :z, :v, :w]
(Multivariate polynomial ring in 5 variables over QQ, QQMPolyRingElem[x, y, z, v, w])

julia> g = QQMPolyRingElem[z^3-z*w^8-z - v^2, 3*z^2-w^8-1, v, z*w, x^8+1, -x^6+x^4-x^2-y^2]

6-element Vector{QQMPolyRingElem}:
 z^3 - z*w^8 - z - v^2
 3*z^2 - w^8 - 1
 v
 z*w
 x^8 + 1
 -x^6 + x^4 - x^2 - y^2

julia> @time radical(ideal(R, g))
 36.195259 seconds (1.77 k allocations: 52.578 KiB)
Ideal generated by
  v
  z
  x^2*y^2 + x^2 + y^2 - 1
  2*x^4 - 4*x^2 + y^4 - 2*y^2 + 3
  4*x^2 + y^6 - y^4 + 7*y^2 - 3
  w^8 + 1

julia> @time radical(ideal(R, g))
 35.533215 seconds (1.77 k allocations: 52.578 KiB)
Ideal generated by
  v
  z
  x^2*y^2 + x^2 + y^2 - 1
  2*x^4 - 4*x^2 + y^4 - 2*y^2 + 3
  4*x^2 + y^6 - y^4 + 7*y^2 - 3
  w^8 + 1

julia> @time radical(ideal(R, g))
 35.559862 seconds (1.77 k allocations: 52.578 KiB)
Ideal generated by
  v
  z
  x^2*y^2 + x^2 + y^2 - 1
  2*x^4 - 4*x^2 + y^4 - 2*y^2 + 3
  4*x^2 + y^6 - y^4 + 7*y^2 - 3
  w^8 + 1

Bot sets of timings come from "hot" sessions, i.e. I'd run the computation before in the session to everything was already compiled.

fingolfin commented 4 months ago

@hannes14 reported that the issue is resolved on his machine 0.23.2. But as noted above (and just re-confirmed) on e.g. tutulla it still takes 40 seconds, with at time 50 processes running, while on my macOS laptop it finishes in 3-4 seconds. That suggests to me that there is still an issue on linux.

Out of curiosity: @HechtiDerLachs is there any change for you?

HechtiDerLachs commented 4 months ago

After updating Oscar for roughly 45 minutes the original example now finishes within roughly 12 seconds for me on Ubuntu. But that is still too much, given that within Singular the computation takes < 1s. It's progress, however. Thanks!

fingolfin commented 2 weeks ago

On Tutualla (a Linux machine with 1.5 TB ram, and two AMD EPYC 9554 64-core CPUs), the minimized example

R, (x,y,z,v,w) = QQ[:x, :y, :z, :v, :w];
g = [z^3-z*w^8-z - v^2, 3*z^2-w^8-1, v, z*w, x^8+1, -x^6+x^4-x^2-y^2];
@time radical(ideal(R, g))

still takes 30 seconds, while on my MacBook it is just 8 seconds -- though note how in May it was only 3.5 seconds (so it got worse).

So I really don't think this issue is fully resolved, esp. since @HechtiDerLachs reports that doing this in Singular directly takes less than a second.

How can we proceed?

JohnAAbbott commented 2 weeks ago

A quick check with help from @hannes14 shows that the slow part is Singular.LibPrimdec.radical. Hans will summon his courage and investigate... sooner or later.