Hi, On 28 May 2011 15:30, SherjilOzair <[email protected]> wrote:
> Hello everyone, > I've been successful in writing the symbolic cholesky decomposition > for sparse matrices in O(n * c**2) time. This is a reasonable order > for sparse systems, but still the performance is not very good. Using > python bulitins, It factors a 100 * 100 Matrix with sparsity 0.57 in > 961 milli-seconds. Using Sympy's numbers, it takes forever or is pain- > stakingly slow for matrices larger than 20 * 20. > In [1] you will find a very simple comparison of Integer, int and mpz. This applies to the rational case as well, just the difference is even bigger. > > I understand why we must integrate groundtypes in matrices to make it > usable. But I don't know how exactly to do it. > > I see that the Matrix constructor currently employs sympify, so it > changes regular ints to Sympy's Integer. I had removed this when I > wanted to test for the python builtins in my DOKMatrix implementation. > > Here's an idea that we can build on. Add a kwarg argument in the > Matrix constructor called dtype, which could takes values like 'gmpy', > 'python', 'sympy', etc. or more detailed, 'int', 'expr', 'poly' etc.. > So that, before putting the element in the matrix, we convert it to > the required dtype. eg. val = gmpy.mpz(val) > > Is it as simple as this, or am I missing something ? > Following sympy.polys design means that you have to employ static typing (all coefficients in a matrix are of the same type, governed by a domain that understands properties of the type). Suppose we have a matrix M over ZZ, then M[0,0] += 1 is well defined and is fast because it requires only one call to domain.convert() (which will exit almost immediately, depending whether ZZ.dtype is Integer, int, mpz or something else). That was simple, but what about M[0,0] += S(1)/2? 1/2 not in ZZ so += may either fail because there is no way to coerce 1/2 to an integer, but it may also figure out a domain for 1/2 (QQ), upgrade the domain in M and proceed. In polys both scenarios can happen depending whether you use DMP (or any other low-level polynomial representation) or low-level APIs of Poly or high-level APIs of Poly (low-level uses the former and high-level uses the later). The main concern in this case is speed (and type checking but it isn't very strong). Deciding whether 1 is in ZZ is fast, but figuring out a domain for 1/2, unifying domain of 1/2 with ZZ (a sup domain has to be found, which in this case is simple, but may be highly non-trivial in case of composite or algebraic domains) and coercing all elements of M, is slow. How to figure out a domain for a set of coefficients? Use construct_domain(). It will give you the domain and will coerce all inputs. Refer to Poly.__new__ and all Poly.from_* and Poly._from_* to see how this works. sympy.polys should have all tools you will need, so try not to reinvent things that are already in SymPy. For example speaking about those "detailed types" 'int', 'expr', 'poly': poly -> what domain and variables?, expr -> what simplification algorithm?, etc. Learn to use what the library provides to you. If there is something missing, e.g. you would like construct_domain() to work with nested lists, that can be done, either on your own or just ask it. For now it may be a little tedious to use stuff from sympy.polys in matrices and at some point I will have share, e.g. domains, with other modules. My suggestion is to start from something simple. You can create a new matrix class that will support the bare minimum of operations to replace Matrix in solve_linear_system(). This new matrix class would support domain construction and type conversions using mechanisms from sympy.polys. Change solve_linear_system() to not use simplify() but rely on the ground types to do the job (solving zero equivalence problem). If this works and is fast, then you can build on top of this. > I would like Mateusz especially to comment on this, and also guide me > and help me learn about the Polys structure. > [1] http://mattpap.github.com/masters-thesis/html/src/internals.html#benchmarking-ground-types > > Regards, > Sherjil Ozair > > -- > You received this message because you are subscribed to the Google Groups > "sympy" group. > To post to this group, send email to [email protected]. > To unsubscribe from this group, send email to > [email protected]. > For more options, visit this group at > http://groups.google.com/group/sympy?hl=en. > > Mateusz -- You received this message because you are subscribed to the Google Groups "sympy" group. To post to this group, send email to [email protected]. To unsubscribe from this group, send email to [email protected]. For more options, visit this group at http://groups.google.com/group/sympy?hl=en.
