GSoC 2014 Application Zamrath Nizam: Efficient Groebner bases and their applications - sympy/sympy GitHub Wiki

Efficient Groebner bases and their applications

Name and Contact Info

Name: Zamrath Nizam

University: University of Moratuwa, Sri Lanka

e-mail: [email protected]

IRC handle: zami

gitter: Zamrath

Skype: zamrath.nizam

Phone No: Will be provided upon request

Abstract

SymPy is an emerging full featured Computer Algebra System (CAS) which provides a platform to deal with symbolic mathematics. This is completely written and used in Python while the code is kept as comprehensible and easily extensible as possible. Groebner basis is a set of non-linear polynomials which provides simple arithmetic solutions for many fundamental problems in the fields of mathematics, natural and technical sciences. In SymPy, this computation uses Buchberger and F5B algorithms but not the most efficient algorithm to date. This proposal is to improve the efficiency of the Groebner base computation using better selection criteria for critical pairs, implement F4 algorithm along with other algorithms and analyze which approach is better in differing contexts, implementing Groebner walk and apply Groebner basis in integration and rational expression simplification.

Community at Large

SymPy is a platform where users can easily perform symbolic mathematical operations such as symbolic integration, differentiation, etc even for multivariate polynomials. Python itself a handy language which comes with rich inbuilt libraries. SymPy comes as another library for Python and it can be run in anywhere a python script can run. Further the comprehensible and easily extensible nature of SymPy provides the open source community a better value on extending their hands.

Groebner basis is one of the most important tool in Linear algebra which is being used in closed form symbolic integration and summation, homogenization/finding projective closure, testing for ideal membership, finding the initial ideal(which gives us lots of information), multivariate polynomial LCM and GCD computation, rational expression manipulation, solving system of equations, etc. Intriguingly this basis solves problems that are far away from algebra as construction of graph colorings, the solution of Sudoku games can be reduced to Groebner bases computation. Reverse Engineering of Software, intelligent control of oil platform, finding genetic relationships between species, solving the inverse robotic kinematics are some example of natural and technical science application [1].

A large variety of problems arising with arbitrary set of polynomials can be solved easily by using equivalent Groebner basis which would otherwise be die-hard. Groebner bases computation is implemented in two algorithms so far namely buchberger and F5B. But these algorithms take considerable amount of resources when the dimensionality of the context (matrices) gets charged up. If, Arbitrary set of polynomial := B, Equivalent Groebner basis := G, Algorithm := buchberger, Total degree of polynomials in B and number of variables:= d, n then the G can be much larger than the B and the upper bound of degree of polynomials in G can be given as 2(1/2*d^2 + d)^(2^(n-1)) which provides an insight on the memory requirement and run-time of the above algorithms.

As par to the original algorithm, the S-polynomial of every pairs of polynomial in the basis should be calculated. To select the pairs intelligently, bruteforce gives nothing at all but Buchberger came up with two elimination criterion named product and chain but later investigation reveled that it can be further brings down while the algorithm is on the process too. These naive algorithms give the way of eliminating one S-polynomial at a time but later Faugere proposed a simultaneous elimination of S-polynomial which accelerated the speed even further under the name F4 [2]. Even further this simultaneous reduction can be transformed in to computation of reduced echelon form problem. This gives advantage of a sparse matrix and techniques in numerical analysis will provide us the efficient way of handling (F5). Solving the min-rank problem for five 6x6 matrices and rank deficiency at least 2, Groebner basis calculation would yield 53 seconds in Buchberger and 1.3 seconds in F5. With 7x7 matrices yield 1000+ seconds and 6.4 seconds respectively[4].

Deliverables and Their Implementation

Selection Criteria

It should be noted that during the execution of Buchberger algorithm, selecting the critical pair and the reductor from a list of polynomials is crucial. It was proven that these choices make no difference in correctness. But the above choices lead to an huge difference in run-time of the computation. This process yields a quick computation as a predefined degree (say d) on which the Buchberger algorithm should be truncated is defined. Here a degree named "Sugar" (say Sd) is calculated for each polynomial. In this sense we reject all the critical pairs whose Sd > d. In practice this algorithm is less efficient. (However the sugar function has been already implemented in sdm_groebner function in distributedmodules.py in polys.)

One thing that should be kept in mind is new algorithms such as F4 has outperformed the Buchberger with any of the selection strategies which have been proposed. Since the inherent efficiency shortfalls of sugar, new versions of Computer Algebra Systems use the Normal strategy explicitly. This strategy is suggested in section 2.5 in [3] where only the critical pairs with a minimal total degree will be accepted. The pseudo code for this is as follows.

sel_normal():
    Input : P a list of critical pairs
    Output: a list of critical pairs
        d  := min{deg(lcm(p)), for all p of P}
        P_d:= {all p of P | deg(lcm(p)) = d}
        return P_d

In groebnertools.py, the selection of critical pairs are based on select function where it takes the pair with minimum LCM(LM(f), LM(g)). But this is less efficient in some contexts such as when all the input polynomials have the same leading term and hence all critical pairs are equal. Above sel_normal can replace the select function since it's proven efficiency over sugar.

F4 Implementation with Analysis on Approaches For Differing Contexts

In contrast to Buchberger, the F4 has more powerful criterion to remove useless critical pairs (S-polynomials). It should be noted that this algorithm does not provide an improvement in worst case scenario. But it is super fast in integer and modulo p polynomial than what can expect in previous implementation. Another factor which rises the speed is; even this algorithm works for any admissible ordering; the DRL ordering which gives the most efficient look for this algorithm. The 'sugar method' implementation on Buchberger gives rise to strong sequential ordering and hence results in difficulty in parallelization. F4 made it possible and provided a way of reducing simultaneously many polynomials utilizing linear algebra techniques. This ends with an analysis in which approach is better in what contexts. The pseudo code for F4[3] is given below.

F4():
    Input : F a finite subset of R[x], Sel a function there exist a Sel(u) for non trivial u
    Output: A finite subset of R[x]

    G, F(+), d := F, F, 0
    P := {Pair(f,g) | f,g of G with f != g}

    while P is not zero set do:
        d   := d + 1
        P_d := Sel(P)
        P   := P\P_d
        L_d := Left(P_d) U Right(P_d)
        F(+):= reduction(L_d,G)

        for h in F(+) do:
           P := P U {Pair(h,g) | for g of G}
           G := G U {h}

In F4, the reduction of a polynomial modulo a subset of R[x] should be extended to the reduction of a subset of R[x] modulo another subset of R[x]. This can be done by calling reduction function which accepts L a finite subset of T x R[x] and G a finite subset of R[x]. The red_groebner function already implemented in groebnertools.py should be extended as illustrated above on behalf of the reduction function in F4. In order to output (for F(+)) the {f of F(hat) whereas HT(f) is not part of HT(F)} where F(hat) is reduction to raw echelon form of F, the symbolic_preprocessing of L and G should be calculated. This is of course the F which is the important factor in this algorithm. The output of symbolic_preprocessing which is F will be a matrix. The symbolic_preprocessing function accepts L and G which are the same inputs provided to the reduction function and output being a finite subset of R[x]. The pseudo code for this is given as follows.

symbolic_preprocessing():

    F := {t*f | (t,f) of L}
    Done := HT(F)

    while T(F) is not Done do:
        m an element of T(F)\Done
        Done := Done U {m}
        if m top reducible modulo G then
            m := m' * HT(f) for some f of G and some m' of T
            F := F U {m' * f}
    return F

Groebner Walk

In some application, a particular ordering of power products would make the Groebner basis computation efficient. In this context the computation should perform in this ordering and transform this Groebner basis to another Groebner basis which is based on the desired ordering. The Groebner walk is a technique which transforms such orderings from one monomial form to another. For ideals of zero dimensional (finite number of solution), fglm algorithm would be enough while for other generic ideals a more extensive algorithm like Groebner walk is needed.

In general if the two regions for the orderings do not overlap, it is always possible to construct a path where all the consecutive orderings in the path have overlaps in their regions. This idea can be better pictured by following example.

Consider the two orderings as <_1 and <_2 where the starting order is the former whereas targeting order is the later. Consider as the Groebner basis is computed with respect to I and <_1 and the C(<_1) and C(<_2) are corresponding Groebner cones. And consider two points s and e as starting point in C(<_1) and ending point in C(<_2) respectively and walk along with the segment. As walk reaches a face of the cone it next enters in to the neighboring cone (say <_new) and first transform the Groebner basis over C(<_1) into the Groebner basis over C(<_new). The fact which ease our life here is the efficient conversion techniques which needs no Buchberger implementation to the Groebner basis over C(<_1). If the reduced Groebner basis over C(<_new) is equivalent to the C(<_2) then the computed Groebner basis over C(<_new) will be the targeted Groebner basis of Ideal I in the ordering of <_2. Otherwise follow the segment until to reach the <_2. This will be terminated in finite number of conversions. The psuedo code for Groebner walk is provided below.

walk():
    Input : G the reduced Groebner basis for I the ideal and (<_1) the ordering
            (<_2) the targeting terms ordering
    Output: G' the reduced Groebner basis for (<_2) the ordering with respect to I the ideal

        w <= -oo
        while compute_last_w() is not oo
            compute generators in_w(G) for in_w(I)
            compute reduced Groebner basis H for in_w(I) over (<_2)
            H' <= {f-f^G | for all f in H}
            Autoreduce H' and put G := H'

        Output <= G

The walk function will be called for the Groebner basis conversion in terms of ordering. To store and compare the vectors of term ordering, a utility class named CallMonomial will be created. This will be able to create instance of specific monomial orderings such as LexOrder, GradedLexOrder, ReversedGradedLexOrder, corresponding InverseOrders and ProductOrder. The walk call will be instantiated with arguments A Groebner basis (G) correspondent to an ideal (I) and the targeting terms ordering ((<_2)). The termination condition for the loop is given by w which is computed by a special function compute_last_w where if w is defined as oo (which represents a considerably large vector set), it will output the G as G'. In compute_last_w, the output will be min_<{of all v of V} where V := {v of (D(G) and C_{<_1,<_2}) for all w<v} and oo if V is empty. As long as w is valid, the generators, reduced Groebner basis H and H' will be computed. Finally the auto_reduce function will be called upon. The auto_reduce functions already implemented in Buchberger will be utilized here.

Applications

This is a work around in some existing code blocks which needs Groebner basis to support their operations such as simplification of rational expression modulo a polynomial ideal and integration functions including rational and transcendental. Moreover continue to extend the functions which needs Groebner bases to support their operations.

Time Schedule

A weekly basis time distribution for my activities is given below. Before the coding phase kick off I have planned to continue the coding and provide a better working source code for this patch. At the mean time I will try to have a more closer insight on already implemented algorithms related to Groebner basis.

May 19 - May 25: Implement the "Normal strategy for F4" to improve the selection criteria
May 26 - Jun 01: Apply this strategy where the old select function have been used. Compare the outputs of undergoing functions and see whether those will give the desired outputs
Jun 02 - Jun 08: Implement F4 algorithm
Jun 09 - Jun 15: Finish and test the source code for F4
Jun 16 - Jun 22: Analyze which approach is better in what contexts. Implement an interface by defining the default Groebner basis computation algorithm as F4 and providing user specification options to switch in to other algorithms as well
Jun 23 - Jun 29: Measure the performance of F4 along with other Groebner basis algorithms and carry out the optimizations
Jun 30 - Jul 06: Implement Groebner walk function with defining utility classes which support functions of walk
Jul 07 - Jul 13: Continue works on Groebner walk
Jul 14 - Jul 20: Implement auto_order function to compromise the term ordering on which Groebner basis should be computed using run-time of computations over every possible ordering with a benchmark
Jul 21 - Jul 27: Define order_options function which carries options to carry on the Groebner basis computation in either a specific term ordering or a heuristic way
Jul 28 - Aug 03: Integrate Groebner walk in to the order_options and check the performance. Apply the Groebner basis in functions which I have mentioned above
Aug 04 - Aug 10: Continue extending applications where the Groebner basis would be highly demanded to support the Groebner basis improvements

At the meantime I would be able to push the changes as soon as I am done with a considerable amount of code as this would make the discussion on working progress more comprehensible. Along with these a document describing functions also will be delivered.

Post GSoC

There are lots of researches underway in this particular domain. Several efficient selection criterion also being looked and tested. After successfully contributed to this project, I have determined to work on merging new proven efficient algorithms over old ones whenever possible.

Me and My Programming Background

I am a third year undergraduate student who is pursuing Electronic and Telecommunication Engineering BSc(Hons) degree in University of Moratuwa, Sri Lanka. This is my first ever attempt to any GSoC version. I have completed modules related to Linear Algebra, Programming, Modular Software Development, Data Structures and Algorithm, etc in the same university mentioned above. Familiar with C, C++, Java, JavaScript and especially with Python. Moreover I have coded some projects in Verilog, MATLAB, etc. Related to above mentioned modules, I could compile a GUI interface for student enrollment system, Library management system, etc. Truth to be told, I am still a beginner level user in version control platforms. Seniors guided me in this and also I went through resources in the web and could get the basic knowledge. Apart from these, have experience in real time programming competitions such as IEEEXtreme, Google CodeJam, FB HackerCup, etc.

Briefly, I am particularly keen on programming stuffs as well as mathematics and hence it would be an intriguing one off experience to deal with these both ends in one platform. I have determined to devote 40 hours per week as to complete this project timely.

Contribution to SymPy

The pull requests I have made related to SymPy are as follows.