Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Improve speed of saturation #4485

Open
wants to merge 9 commits into
base: master
Choose a base branch
from

Conversation

paemurru
Copy link
Contributor

@paemurru paemurru commented Jan 18, 2025

The function remove seems to be faster than Singular.saturation.

Benchmark 1

julia> function random_polynomial(
           vars::AbstractVector{T},
           deg::Int64,
           coeff_range::AbstractVector
       ) where T<:MPolyDecRingElem
         S = parent(vars[1])
         indices = Int64[]
         for var in vars
           append!(indices, [findfirst(isequal(var), gens(S))])
         end
         monomials = collect(monomials_of_degree(S, deg, indices))
         coefficients = rand(coeff_range, length(monomials))
         return sum(map(*, monomials, coefficients))
       end
random_polynomial (generic function with 2 methods)

julia> R, x = graded_polynomial_ring(QQ, "x" => 1:10)
(Graded multivariate polynomial ring in 10 variables over QQ, MPolyDecRingElem{QQFieldElem, QQMPolyRingElem}[x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10]])

julia> f = random_polynomial(gens(R), 10, 0:1);

Before:

julia> @btime saturation(ideal(f), ideal(x[1]));
  3.585 s (3780311 allocations: 79.62 MiB)

Now:

julia> @btime saturation(ideal(f), ideal(x[1]));
  10.827 ms (92416 allocations: 9.35 MiB)

Benchmark 2

The line phi_star = Oscar.pushforward_on_algebraic_lattices(phi) in the testset "moebius transformations" in the tests of AlgebraicGeometry/Schemes/EllipticSurface.

Before:

julia> @btime phi_star = Oscar.pushforward_on_algebraic_lattices(phi)
  20.848 s (54676949 allocations: 2.21 GiB)
Map
  from quadratic space of dimension 10
  to quadratic space of dimension 10

Now:

julia> @btime phi_star = Oscar.pushforward_on_algebraic_lattices(phi)
  2.800 s (17432516 allocations: 609.05 MiB)
Map
  from quadratic space of dimension 10
  to quadratic space of dimension 10

@ederc
Copy link
Member

ederc commented Jan 20, 2025

For one variable I can see that this helps keeping the overhead of calling Singular's saturation down, but I suspect that for saturating w.r.t. a general polynomial you might be able to construct examples, where remove is not faster.

@paemurru
Copy link
Contributor Author

saturating w.r.t. a general polynomial you might be able to construct examples, where remove is not faster.

Before the last commit, some tests did not finish (https://github.com/oscar-system/Oscar.jl/actions/runs/12854953830), so I think you are right

@simonbrandhorst
Copy link
Collaborator

simonbrandhorst commented Jan 20, 2025

Playing a bit with random examples for f and g. I get consistently a factor of at least 100x speedup of remove over saturation for principal ideals. Do not think this is just the call overhead of singular saturation. It seems impossible to me that the call overhead is 30 seconds.
In some examples singular will not even terminate but remove terminates in a fraction of a second.

julia> f = random_polynomial(gens(R)[1:5], 10, 0:100);

julia> g = random_polynomial(gens(R)[1:5], 10, 0:1);

julia> @time remove(J[1],g);
  0.009453 seconds (19 allocations: 45.336 KiB)

julia> @time saturation(J,ideal(g));
 32.858295 seconds (121.48 M allocations: 26.022 GiB, 34.39% gc time)

julia> f = random_polynomial(gens(R)[1:5], 10, 0:100);

julia> g = random_polynomial(gens(R)[1:5], 10, 0:1);

julia> J = ideal(g*f);

julia> @time remove(J[1],g);
  0.009289 seconds (19 allocations: 45.078 KiB)

julia> @time saturation(J,ideal(g));
  1.988479 seconds (9.95 M allocations: 750.730 MiB, 32.38% gc time)

Does singular catch the special case of principal ideals in its preprocessing?
My impression is that a generic algorithm is used.

Can anyone produce an example where singular is faster?

@fingolfin fingolfin requested review from wdecker and ederc and removed request for wdecker January 20, 2025 22:10
src/Rings/mpoly-ideals.jl Outdated Show resolved Hide resolved
src/Rings/mpoly-ideals.jl Outdated Show resolved Hide resolved
Copy link
Member

@ederc ederc left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

After running a couple of tests I think using remove for the principal ideal case is really worthwhile, not only if J = <x> for a generator x.

Still we need to fix the tests.

@ederc
Copy link
Member

ederc commented Jan 21, 2025

The failing book test is just a sign difference:

<     2: scheme((z//x)^3*t^7 - 2*(z//x)^3*t^6 + (z//x)^3*t^5 - (z//x)*(y//x)^2 + 1)
<     3: scheme((z//y)^3*t^7 - 2*(z//y)^3*t^6 + (z//y)^3*t^5 - (z//y) + (x//y)^3)
<     4: scheme((x//z)^3 - (y//z)^2 + s^7 - 2*s^6 + s^5)
<     5: scheme((z//x)^3*s^7 - 2*(z//x)^3*s^6 + (z//x)^3*s^5 - (z//x)*(y//x)^2 + 1)
<     6: scheme((z//y)^3*s^7 - 2*(z//y)^3*s^6 + (z//y)^3*s^5 - (z//y) + (x//y)^3)
---
>     2: scheme(-(z//x)^3*t^7 + 2*(z//x)^3*t^6 - (z//x)^3*t^5 + (z//x)*(y//x)^2 - 1)
>     3: scheme(-(z//y)^3*t^7 + 2*(z//y)^3*t^6 - (z//y)^3*t^5 + (z//y) - (x//y)^3)
>     4: scheme(-(x//z)^3 + (y//z)^2 - s^7 + 2*s^6 - s^5)
>     5: scheme(-(z//x)^3*s^7 + 2*(z//x)^3*s^6 - (z//x)^3*s^5 + (z//x)*(y//x)^2 - 1)
>     6: scheme(-(z//y)^3*s^7 + 2*(z//y)^3*s^6 - (z//y)^3*s^5 + (z//y) - (x//y)^3)

But I am not sure why the other tests take too long and get cancelled.

Copy link
Member

@benlorenz benlorenz left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks like all the long tests are timing out, getting stuck in some computation from test/AlgebraicGeometry/Schemes/EllipticSurface.jl.
Probably the moebius transformations test-set, this should appear after two neighbor step but doesn't appear in the logs.
I have started that test-file locally but it is still running after 20 minutes (while on CI it takes about 5).
Sending USR1 does indeed point to a remove call:

signal (10): User defined signal 1
_nmod_mpoly_divides_monagan_pearce at /workspace/srcdir/flint/src/nmod_mpoly/divides_monagan_pearce.c:232
nmod_mpoly_divides_monagan_pearce at /workspace/srcdir/flint/src/nmod_mpoly/divides_monagan_pearce.c:596
nmod_mpoly_divides at /workspace/srcdir/flint/src/nmod_mpoly/divides.c:192
divides at /home/lorenz/.julia/packages/Nemo/y7UoI/src/flint/nmod_mpoly.jl:483
unknown function (ip: 0x7f89d711167a)
divides at /home/lorenz/.julia/packages/Nemo/y7UoI/src/flint/fq_default_mpoly.jl:339
remove at /home/lorenz/.julia/packages/AbstractAlgebra/UzHEi/src/MPoly.jl:839
#saturation#305 at /home/lorenz/software/polymake/julia/Oscar.jl/src/Rings/mpoly-ideals.jl:352
saturation at /home/lorenz/software/polymake/julia/Oscar.jl/src/Rings/mpoly-ideals.jl:349
unknown function (ip: 0x7f89d69bd4e6)
saturation at /home/lorenz/software/polymake/julia/Oscar.jl/src/Rings/mpolyquo-localizations.jl:2197
_iterative_saturation at /home/lorenz/software/polymake/julia/Oscar.jl/src/AlgebraicGeometry/Schemes/Sheaves/IdealSheaves.jl:571
_iterative_saturation at /home/lorenz/software/polymake/julia/Oscar.jl/src/AlgebraicGeometry/Schemes/Sheaves/IdealSheaves.jl:565
unknown function (ip: 0x7f89c8d0ac96)
#produce_object#5672 at /home/lorenz/software/polymake/julia/Oscar.jl/src/AlgebraicGeometry/Schemes/Sheaves/IdealSheaves.jl:1792
produce_object at /home/lorenz/software/polymake/julia/Oscar.jl/src/AlgebraicGeometry/Schemes/Sheaves/IdealSheaves.jl:1696
#_#5306 at /home/lorenz/software/polymake/julia/Oscar.jl/src/AlgebraicGeometry/Schemes/Sheaves/Sheaves.jl:52
AbsPreSheaf at /home/lorenz/software/polymake/julia/Oscar.jl/src/AlgebraicGeometry/Schemes/Sheaves/Sheaves.jl:50 [inlined]
#5691 at ./none:0
unknown function (ip: 0x7f89c8f5a8d2)
iterate at ./generator.jl:48 [inlined]
collect_to! at ./array.jl:849
collect_to_with_first! at ./array.jl:827
unknown function (ip: 0x7f89c8f598e5)
collect at ./array.jl:801
unknown function (ip: 0x7f89c8f59392)
produce_object at /home/lorenz/software/polymake/julia/Oscar.jl/src/AlgebraicGeometry/Schemes/Sheaves/IdealSheaves.jl:1819
#_#5306 at /home/lorenz/software/polymake/julia/Oscar.jl/src/AlgebraicGeometry/Schemes/Sheaves/Sheaves.jl:52
AbsPreSheaf at /home/lorenz/software/polymake/julia/Oscar.jl/src/AlgebraicGeometry/Schemes/Sheaves/Sheaves.jl:50
unknown function (ip: 0x7f89c8f58ff2)
#colength#5764 at /home/lorenz/software/polymake/julia/Oscar.jl/src/AlgebraicGeometry/Schemes/Sheaves/IdealSheaves.jl:2288
colength at /home/lorenz/software/polymake/julia/Oscar.jl/src/AlgebraicGeometry/Schemes/Sheaves/IdealSheaves.jl:2259
unknown function (ip: 0x7f89c8d04cc6)
#intersect#5930 at /home/lorenz/software/polymake/julia/Oscar.jl/src/AlgebraicGeometry/Schemes/Divisors/WeilDivisor.jl:344
intersect at /home/lorenz/software/polymake/julia/Oscar.jl/src/AlgebraicGeometry/Schemes/Divisors/WeilDivisor.jl:320
unknown function (ip: 0x7f89c2aa3ad6)
basis_representation at /home/lorenz/software/polymake/julia/Oscar.jl/src/AlgebraicGeometry/Surfaces/EllipticSurface/EllipticSurface.jl:994
#7050 at ./none:0
unknown function (ip: 0x7f89c2aa3772)
iterate at ./generator.jl:48 [inlined]
collect at ./array.jl:791
unknown function (ip: 0x7f89c2aa13a2)
pushforward_on_algebraic_lattices at /home/lorenz/software/polymake/julia/Oscar.jl/src/AlgebraicGeometry/Surfaces/EllipticSurface/Morphisms.jl:430
unknown function (ip: 0x7f89c2a2b2c2)
jl_apply at /cache/build/tester-amdci5-12/julialang/julia-release-1-dot-11/src/julia.h:2157 [inlined]
do_call at /cache/build/tester-amdci5-12/julialang/julia-release-1-dot-11/src/interpreter.c:126
eval_value at /cache/build/tester-amdci5-12/julialang/julia-release-1-dot-11/src/interpreter.c:223
eval_body at /cache/build/tester-amdci5-12/julialang/julia-release-1-dot-11/src/interpreter.c:562
eval_body at /cache/build/tester-amdci5-12/julialang/julia-release-1-dot-11/src/interpreter.c:539
eval_body at /cache/build/tester-amdci5-12/julialang/julia-release-1-dot-11/src/interpreter.c:539
jl_interpret_toplevel_thunk at /cache/build/tester-amdci5-12/julialang/julia-release-1-dot-11/src/interpreter.c:821
jl_toplevel_eval_flex at /cache/build/tester-amdci5-12/julialang/julia-release-1-dot-11/src/toplevel.c:943
jl_toplevel_eval_flex at /cache/build/tester-amdci5-12/julialang/julia-release-1-dot-11/src/toplevel.c:886
ijl_toplevel_eval_in at /cache/build/tester-amdci5-12/julialang/julia-release-1-dot-11/src/toplevel.c:994
eval at ./boot.jl:430 [inlined]
include_string at ./loading.jl:2734
_include at ./loading.jl:2794
include at ./Base.jl:558 [inlined]
macro expansion at ./timing.jl:581 [inlined]
#_timed_include#33 at /home/lorenz/software/polymake/julia/Oscar.jl/src/utils/tests.jl:5
_timed_include at /home/lorenz/software/polymake/julia/Oscar.jl/src/utils/tests.jl:1 [inlined]
_timed_include at /home/lorenz/software/polymake/julia/Oscar.jl/src/utils/tests.jl:1 [inlined]
#40 at /home/lorenz/software/polymake/julia/Oscar.jl/src/utils/tests.jl:258
with_unicode at /home/lorenz/.julia/packages/AbstractAlgebra/UzHEi/src/PrettyPrinting.jl:1846
#test_module#39 at /home/lorenz/software/polymake/julia/Oscar.jl/src/utils/tests.jl:206 [inlined]
test_module at /home/lorenz/software/polymake/julia/Oscar.jl/src/utils/tests.jl:205
unknown function (ip: 0x7f8a1538af66)
jl_apply at /cache/build/tester-amdci5-12/julialang/julia-release-1-dot-11/src/julia.h:2157 [inlined]
do_call at /cache/build/tester-amdci5-12/julialang/julia-release-1-dot-11/src/interpreter.c:126
eval_value at /cache/build/tester-amdci5-12/julialang/julia-release-1-dot-11/src/interpreter.c:223
eval_stmt_value at /cache/build/tester-amdci5-12/julialang/julia-release-1-dot-11/src/interpreter.c:174 [inlined]
eval_body at /cache/build/tester-amdci5-12/julialang/julia-release-1-dot-11/src/interpreter.c:663
jl_interpret_toplevel_thunk at /cache/build/tester-amdci5-12/julialang/julia-release-1-dot-11/src/interpreter.c:821
jl_toplevel_eval_flex at /cache/build/tester-amdci5-12/julialang/julia-release-1-dot-11/src/toplevel.c:943
jl_toplevel_eval_flex at /cache/build/tester-amdci5-12/julialang/julia-release-1-dot-11/src/toplevel.c:886
ijl_toplevel_eval_in at /cache/build/tester-amdci5-12/julialang/julia-release-1-dot-11/src/toplevel.c:994
eval at ./boot.jl:430 [inlined]
eval_user_input at /cache/build/tester-amdci5-12/julialang/julia-release-1-dot-11/usr/share/julia/stdlib/v1.11/REPL/src/REPL.jl:245
repl_backend_loop at /cache/build/tester-amdci5-12/julialang/julia-release-1-dot-11/usr/share/julia/stdlib/v1.11/REPL/src/REPL.jl:342
#start_repl_backend#59 at /cache/build/tester-amdci5-12/julialang/julia-release-1-dot-11/usr/share/julia/stdlib/v1.11/REPL/src/REPL.jl:327
start_repl_backend at /cache/build/tester-amdci5-12/julialang/julia-release-1-dot-11/usr/share/julia/stdlib/v1.11/REPL/src/REPL.jl:324
#run_repl#72 at /cache/build/tester-amdci5-12/julialang/julia-release-1-dot-11/usr/share/julia/stdlib/v1.11/REPL/src/REPL.jl:483
run_repl at /cache/build/tester-amdci5-12/julialang/julia-release-1-dot-11/usr/share/julia/stdlib/v1.11/REPL/src/REPL.jl:469
jfptr_run_repl_10104.1 at /home/lorenz/software/polymake/julia/julia/julia-1.11.2/share/julia/compiled/v1.11/REPL/u0gqU_4x0TT.so (unknown line)
#1150 at ./client.jl:446
jfptr_YY.1150_14803.1 at /home/lorenz/software/polymake/julia/julia/julia-1.11.2/share/julia/compiled/v1.11/REPL/u0gqU_4x0TT.so (unknown line)
jl_apply at /cache/build/tester-amdci5-12/julialang/julia-release-1-dot-11/src/julia.h:2157 [inlined]
jl_f__call_latest at /cache/build/tester-amdci5-12/julialang/julia-release-1-dot-11/src/builtins.c:875
#invokelatest#2 at ./essentials.jl:1055 [inlined]
invokelatest at ./essentials.jl:1052 [inlined]
run_main_repl at ./client.jl:430
repl_main at ./client.jl:567 [inlined]
_start at ./client.jl:541
jfptr__start_73406.1 at /home/lorenz/software/polymake/julia/julia/julia-1.11.2/lib/julia/sys.so (unknown line)
jl_apply at /cache/build/tester-amdci5-12/julialang/julia-release-1-dot-11/src/julia.h:2157 [inlined]
true_main at /cache/build/tester-amdci5-12/julialang/julia-release-1-dot-11/src/jlapi.c:900
jl_repl_entrypoint at /cache/build/tester-amdci5-12/julialang/julia-release-1-dot-11/src/jlapi.c:1059
main at julia-latest (unknown line)
unknown function (ip: 0x7f8a1d3e047f)
__libc_start_main at /lib64/libc.so.6 (unknown line)
unknown function (ip: 0x4010b8)
unknown function (ip: (nil))

So maybe looking into whatever ideals are generated there might produce examples where saturation is faster?

src/Rings/mpoly-ideals.jl Outdated Show resolved Hide resolved
src/Rings/mpoly-ideals.jl Outdated Show resolved Hide resolved
@paemurru
Copy link
Contributor Author

It is the function Oscar.pushforward_on_algebraic_lattices from Benchmark 2 in the top post that does not finish. For easier debugging, I have added the code here

k = GF(113)
kt,t = polynomial_ring(k,:t)
Ft = fraction_field(kt)

E = elliptic_curve(Ft, Ft.([0,0,0,(56*t^8 + 56),0]));
basis_mwl = [[112*t^4 + 112*t^3 + 56*t^2 + 15*t + 1, 100*t^6 + 24*t^5 + 56*t^4 + 13*t^3 + 64*t^2 + 89*t + 31],
             [31*t^4 + 15*t^2 + 82, 44*t^5 + 14*t^3 + 69*t],
             [82*t^4 + 13, 37*t^4 + 10],
             [91*t^4 + 16*t^3 + 25*t^2 + 14*t + 22, 18*t^6 + 55*t^5 + 45*t^4 + 44*t^3 + 110*t^2 + 58*t + 69],
             [32*t^4 + 72*t^2 + 81, 78*t^6 + 79*t^4 + 58*t^2 + 40],
             [56*t^4 + 49*t^3 + 85*t^2 + 57*t + 57, 72*t^6 + 25*t^5 + 22*t^4 + 101*t^3 + 104*t^2 + 88*t + 50],
             [22*t^4 + 87*t^3 + 77*t^2 + 26*t + 22, 44*t^6 + 27*t^5 + 68*t^4 + 98*t^3 + 45*t^2 + 27*t + 69],
             [112*t^4 + 44*t^3 + 49*t^2 + 44*t + 112, 100*t^6 + 74*t^5 + 49*t^4 + 8*t^3 + 49*t^2 + 74*t + 100]]
basis_mwl = [E(i) for i in basis_mwl];  

X = elliptic_surface(E, 2, basis_mwl[1:0])  # speed up test by hiding the sections

moeb = Oscar.admissible_moebius_transformations(X, X)

@test length(moeb) == 16 # This must really be the number for this particular surface.

for phi in moeb
  @test Oscar.check_isomorphism_on_generic_fibers(phi)
end

phi = moeb[5] # Just pick one which is not the identity

W = weierstrass_chart_on_minimal_model(X)

for (i, U) in enumerate(affine_charts(X))
  phi_loc = Oscar.cheap_realization(phi, W, U)
  @test all(iszero(pullback(phi_loc)(lifted_numerator(g))) for g in gens(modulus(OO(U))))
end

phi_star = Oscar.pushforward_on_algebraic_lattices(phi)

@ederc
Copy link
Member

ederc commented Jan 22, 2025

Here is a minimal example of the problem:

R = @polynomial_ring(QQ, [:x])
I = ideal(R,1)
J = ideal(R,1)
remove(I[1],J[1])

remove runs into an endless loop calling divides on 1 and 1.
So either we adjust remove in AbstractAlgebra to allow such input or we catch this special case in Oscar.

@thofma
Copy link
Collaborator

thofma commented Jan 22, 2025

Here is a minimal example of the problem:

R = @polynomial_ring(QQ, [:x])
I = ideal(R,1)
J = ideal(R,1)
remove(I[1],J[1])

remove runs into an endless loop calling divides on 1 and 1. So either we adjust remove in AbstractAlgebra to allow such input or we catch this special case in Oscar.

For univariates we have

julia> remove(x, x^0)
ERROR: Second argument must be a non-zero non-unit
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35

I will adjust it in AA to throw the same error for the multivariates. The usage here should be adjusted.

Edit: Nemocas/AbstractAlgebra.jl#1967

@paemurru
Copy link
Contributor Author

Adjusted the code here like @thofma suggested. Assuming the tests pass, I believe this can now be merged.

@simonbrandhorst
Copy link
Collaborator

Adjusted the code here like @thofma suggested. Assuming the tests pass, I believe this can now be merged.

Just the signs in the booktests need to be changed.

Copy link

codecov bot commented Jan 23, 2025

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 84.55%. Comparing base (767415a) to head (010dee7).
Report is 10 commits behind head on master.

Additional details and impacted files
@@            Coverage Diff             @@
##           master    #4485      +/-   ##
==========================================
+ Coverage   84.39%   84.55%   +0.16%     
==========================================
  Files         672      672              
  Lines       88690    88887     +197     
==========================================
+ Hits        74846    75156     +310     
+ Misses      13844    13731     -113     
Files with missing lines Coverage Δ
src/Rings/mpoly-ideals.jl 94.22% <100.00%> (+0.07%) ⬆️

... and 25 files with indirect coverage changes

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

6 participants