One Hat Cyber Team
Your IP :
216.73.216.115
Server IP :
194.44.31.54
Server :
Linux zen.imath.kiev.ua 4.18.0-553.77.1.el8_10.x86_64 #1 SMP Fri Oct 3 14:30:23 UTC 2025 x86_64
Server Software :
Apache/2.4.37 (Rocky Linux) OpenSSL/1.1.1k
PHP Version :
5.6.40
Buat File
|
Buat Folder
Eksekusi
Dir :
~
/
usr
/
share
/
Macaulay2
/
Edit File:
EigenSolver.m2
newPackage( "EigenSolver", Version => "0.1", Date => "June 2020", Authors => { {Name => "Laurent Busé", Email => "Laurent.Buse@inria.fr"}, {Name => "Justin Chen", Email => "justin.chen@math.gatech.edu"}, {Name => "Kisun Lee", Email => "kil004@ucsd.edu"}, {Name => "Anton Leykin", Email => "anton.leykin@gmail.com"}, {Name => "Tomas Pajdla", Email => "pajdla@cvut.cz"}, {Name => "Erika Pirnes", Email => "pirnes@wisc.edu"} }, Headline => "polynomial system solver via eigen-computations", Keywords => {"Numerical Algebraic Geometry"}, PackageExports => {"NAGtypes"}, AuxiliaryFiles => false, DebuggingMode => false ) export { "zeroDimSolve", "Basis", "Multiplier" } zeroDimSolve = method(Options => { symbol Basis => null, symbol Multiplier => 0, Strategy => "Stickelberger", Tolerance => 0.000001 }) zeroDimSolve Ideal := List => opts -> I -> ( -- Assume I is radical if opts.Strategy == "Stickelberger" then return eigSolve1(I, opts) else if opts.Strategy == "elimMatrixP1P1" then return eigSolveP1P1(I, opts) else if opts.Strategy == "elimMatrixP2" then return eigSolveP2(I, opts) else error "Unknown strategy"; ) reduceCoeffs = (A, I, m) -> transpose sub(last coefficients(A % I, Monomials => m), coefficientRing ring I) eigSolve1 = method(Options => options zeroDimSolve) -- based on Tomas' function SolvePolyEqByEigV eigSolve1 Ideal := List => opts -> I -> ( R := ring I; B := if opts.Basis =!= null then opts.Basis else sub(basis(R/I), R); f := if opts.Multiplier =!= 0 then opts.Multiplier else random(1, R); M := reduceCoeffs(f*B, I, B); (E, V) := eigenvectors M; V = V*inverse diagonalMatrix V^{0}; (entries transpose(reduceCoeffs(vars R, I, B)*V))/(p -> point{p}) ) eigSolveP1P1 = method(Options => options zeroDimSolve) -- based on Laurent's code eigSolveP1P1 Ideal := List => opts -> J -> ( -- currently assumes dim R = 2 et 2 equations, homogenizes to P^1 x P^1 R := ring J; FF := coefficientRing R; Lvar := vars R; (xx0, yy0) := (symbol xx0, symbol yy0); nR := FF[xx0,Lvar_(0,0),yy0,Lvar_(0,1)]; E0 := homogenize(homogenize(sub(J_0,nR),xx0,{1,1,0,0}),yy0,{0,0,1,1}); E1 := homogenize(homogenize(sub(J_1,nR),xx0,{1,1,0,0}),yy0,{0,0,1,1}); JnR := ideal(E0,E1); -- homogenized equations nnR := FF[xx0,Lvar_(0,0),yy0,Lvar_(0,1),Degrees=>{{1,0},{1,0},{0,1},{0,1}}]; JnnR := sub(JnR,nnR); -- homogenized equations in the bigraded ring deg := degrees JnnR; satind := {(deg_0)_0+(deg_1)_0-1,(deg_0)_1+(deg_1)_1-1}; psi := map(nnR^{{0,0}}, nnR^(-(JnnR_*/degree)), gens JnnR); elimMat := matrix basis(satind,psi); elimMatTr := transpose sub(elimMat,CC); numrk := numericalRank elimMatTr; K := numericalKernel(elimMatTr,Tolerance=>opts.Tolerance); numroots := rank source K; -- expected number of roots D0 := K^{0..numroots-1}; -- indexed by the basis 1,y_1,.. D1 := K^{satind_1+1..satind_1+numroots}; -- indexed by the basis x_1*1,x_1*y_1,... (EVal,EVec) := eigenvectors (inverse(D0)*D1); -- NEED GENERALIZED EIGENVALUES EVect := D0*EVec; apply(#EVal, i -> point{{EVal_i,EVect_(1,i)/EVect_(0,i)}}) -- roots (x_1,y_1) assuming x0=y0=1 ) eigSolveP2 = method(Options => options zeroDimSolve) -- based on Laurent's code eigSolveP2 Ideal := List => opts -> J -> ( -- currently assumes dim R = 2 et 2 (nonlinear) equations, homogenizes to P^2 R := ring J; FF := coefficientRing R; Lvar := vars R; xx0 := symbol xx0; nR := FF[xx0,gens R]; E0 := homogenize(sub(J_0,nR),xx0); E1 := homogenize(sub(J_1,nR),xx0); JnR := ideal(E0,E1); -- homogenized equations deg := degrees JnR; satind := (deg_0)_0+(deg_1)_0-1; psi := map(nR^{0}, nR^(-(JnR_*/degree)), gens JnR); elimMat := matrix basis(satind,psi); elimMatTr := transpose sub(elimMat,CC); K := numericalKernel(elimMatTr,Tolerance=>opts.Tolerance); numroots := rank source K; -- expected number of roots B := flatten entries basis(satind, nR); H := hashTable apply(#B, i -> B#i => i); preB0 := select(B, b -> member(nR_1*b // nR_0, B)); B0 := {first preB0}; D0 := K^{H#(B0_0)}; -- assumption here i := 0; -- initialization -- to be optimized while (numericalRank(D0) < numroots) and (i < numrows K) do ( i=i+1; b := preB0_i; addRow := K^{H#b}; if numericalRank (D0||addRow) > numericalRank(D0) then ( D0=D0||addRow; B0=B0|{b}; ); ); --B0 := select(numroots, B, b -> member(nR_1*b // nR_0, B)); -- !! must be numerically full rank --D0 := K^(apply(B0, b -> H#b)); B1 := apply(B0, b -> (b*nR_1)//nR_0); D1 := K^(apply(B1, b -> H#b)); (EVal,EVec) := eigenvectors (inverse(D0)*D1); -- NEED GENERALIZED EIGENVALUES (provisionally assume no solutions at infinity, separated by x_1) EVect := D0*EVec; << "numerical rank of D0=" << numericalRank D0 << endl; apply(#EVal, i -> point{{EVal_i,EVect_(2,i)/EVect_(0,i)}}) -- BE CAREFUL: we are assuming that the three rows are indexed by 1,x1,x2. ) -- needs "laurent-eigensolving.m2" -- needs "documentation.m2" TEST /// -- Problem 2.11 in https://math.berkeley.edu/~bernd/cbms.pdf n = 3 R = QQ[a_1..a_(n-1),b_1..b_(n-1)] S = R[x,y]/ideal(x*y-1) f = sum(n-1, i -> R_i*x^(i+1) + R_(n-1+i)*y^(i+1)) + x^n + y^n I = sub(ideal apply(toList(2..2*n-1), i -> last coefficients(f^i, Monomials => {1_S})), R) sols = zeroDimSolve(I, Multiplier => (3/8)*a_1+(3/7)*a_2+(1/3)*b_1+(5/8)*b_2) realPts = realPoints sols assert(#sols == 66 and #realPts == 6) sort apply(sols, p -> norm evaluate(gens I, p)) /// TEST /// needsPackage "ExampleSystems" cyc6 = ideal cyclic(6, QQ) elapsedTime sols = zeroDimSolve cyc6; sort apply(sols, p -> norm evaluate(gens cyc6, p)) /// TEST /// A = QQ[x0,x1,y0,y1,Degrees=>{{1,0},{1,0},{0,1},{0,1}}] I = ideal(random({2,5},A),random({2,5},A)) R = QQ[x1,y1] spe = map(R,A,matrix{{1,x1,1,y1}}) J=spe I listSol = zeroDimSolve(J, Strategy => "elimMatrixP1P1") assert(#listSol == 20) B = QQ[x0,x1,x2] I = ideal(random(2,B), random(2,B)) R = QQ[x1,x2] spe = map(R,B,matrix{{1,x1,x2}}) J = spe I elapsedTime listSol = zeroDimSolve(J, Strategy => "elimMatrixP2") sort apply(listSol, p -> norm evaluate(gens J, p)) zSol= zeroDimSolve(J, Strategy => "Stickelberger") netList zSol sort apply(zSol, p -> norm evaluate(gens J, p)) assert(#listSol == product(I_*/degree/first)) /// TEST /// -- Example from Roser and Olga needsPackage "GraphicalModels" G=graph{{1,2},{2,3},{3,4},{1,4}} R=gaussianRing(G) K=undirectedEdgesMatrix R d=4 numS=lift(d*(d+1)/2,ZZ); lpRvar=apply(numgens(R)-numS,i->(gens(R))_i); R2=coefficientRing(R)[lpRvar] R2map=apply(numgens(R),i-> if i<= numgens(R)-numS-1 then (gens(R2))_i else 0) F=map(R2,R,R2map) K2=F(K) X=random(ZZ^4,ZZ^4) S=X*transpose(X) I=ideal{jacobian ideal{determinant(K2)}-determinant(K2)*jacobian(ideal{trace(K2*S)})} J=saturate(I,ideal{determinant(K2)}) sols = zeroDimSolve J assert(#sols == 5) /// beginDocumentation() doc /// Key EigenSolver Headline polynomial system solver via eigen-computations Description Text This package implements a solver for zero-dimensional polynomial systems based on eigenvalue/eigenvector computations. The theoretical basis for this is Stickelberger's Theorem (cf. [1, Theorem 2.6], also [2]), which states that if I is a zero-dimensional ideal in a polynomial ring $R$ over an algebraically closed field $k$, then the points of V(I) can be obtained as eigenvalues of multiplication matrices corresponding to variables, on the finite-dimensional $k$-vector space $R/I$. Since the main solving is done via linear algebra, this solver can be significantly quicker than other solvers performing nonlinear computations. However, a Grobner basis of the ideal I is still needed, e.g. in order to obtain the sizes of the multiplication matrices. Example needsPackage "ExampleSystems" I = ideal cyclic(6, QQ) elapsedTime sols = zeroDimSolve I; #sols -- 156 solutions Text The authors would like to acknowledge the June 2020 Macaulay2 workshop held virtually at Warwick, where this package was first developed. Text {\bf References}: Code UL { "[1] Sturmfels, Bernd. Solving systems of polynomial equations. No. 97. American Mathematical Soc., 2002", "[2] Cox, David A. Stickelberger and the Eigenvalue Theorem. https://arxiv.org/abs/2007.12573" } /// doc /// Key zeroDimSolve (zeroDimSolve, Ideal) Basis Multiplier Headline zero-dimensional polynomial system solver Usage zeroDimSolve(I) Inputs I:Ideal a zero-dimensional ideal Outputs :List of points on V(I) Description Text This function is a general-purpose solver for zero-dimensional polynomial systems. The default {\tt Strategy} is {\tt "Stickelberger"}, which leverages Stickelberger's theorem to compute the coordinates of solution points as (joint) eigenvalues of multiplication matrices corresponding to the variables of the ring. The option {\tt Multiplier} allows the user to specify the linear form used to construct the multiplication matrix. For example, if the ambient polynomial ring is $QQ[x,y,z]$ and only the $y$-coordinates of solutions are needed, then one can specify {\tt Multiplier => y} (and ignore the coordinates of other solution points). By default, a random linear form is used: due to this, it may be helpful to run this function a couple of times to minimize the likelihood of a "bad" linear form being chosen. The option {\tt Basis} allows the user to specify a basis of the quotient $R/I$, if one is known. The default value is null, in which case a basis for $R/I$ will be computed via a Grobner basis computation. In general, {\tt Strategy => "Stickelberger"} does not require the input ideal I to be radical, so solutions may appear with multiplicity: Example R = QQ[x,y] I = ideal"x2,xy,y3" sols = zeroDimSolve I #sols == 4 and all(sols, p -> clean(1e-16, matrix p) == 0) SeeAlso "solveSystem" /// end-- restart needsPackage "EigenSolver" loadPackage("EigenSolver", Reload => true) installPackage("EigenSolver", RemakeAllDocumentation => true) check "EigenSolver"
Simpan