On a similar note, I haven't seen the RationalFunction type in Sympy which should essentially be Poly/Poly. Poly/Poly current returns a Mul, which would not be very good for the zero-equivalence problem, I think.
Matrix factorization algorithms employ division on elements often. So, if the elements belong to Poly, then we expect the matrix factorization to contain RationalFunction (RationalPoly ?). I think it would be the most natural type for elements of Matrices. It will use the Poly class as its numerator and denominator, thus taking advantage of the well-implemented Poly. Zero-equivalence problem would be to check if the numerator is zero only, etc. What does the community think about this class, the RationalFunction ? Do add ideas from your knowledge and experience also. On May 29, 3:49 pm, SherjilOzair <sherjiloz...@gmail.com> wrote: > Could you explain why normal multiplication/addition operations on > Poly is slower than on Muls and Adds, when Poly is supposed to be > lower level than Mul, Add ? > > In [35]: A = randMatrix(5,5) > > In [36]: A = A.applyfunc(lambda i: i + x) > > In [37]: %timeit B = A.T * A > 1000 loops, best of 3: 1.63 ms per loop > > In [38]: C = A.applyfunc(poly) > > In [39]: C > Out[39]: > ⎡Poly(x + 56, x, domain='ZZ') Poly(x + 89, x, domain='ZZ') Poly(x + > 23, x, domain='ZZ') Poly(x + 46, x, domain='ZZ') Poly(x + 95, x, > domain='ZZ')⎤ > ⎢ > ⎥ > ⎢Poly(x + 41, x, domain='ZZ') Poly(x + 3, x, domain='ZZ') Poly(x + > 17, x, domain='ZZ') Poly(x + 13, x, domain='ZZ') Poly(x + 85, x, > domain='ZZ')⎥ > ⎢ > ⎥ > ⎢Poly(x + 10, x, domain='ZZ') Poly(x + 95, x, domain='ZZ') Poly(x + > 60, x, domain='ZZ') Poly(x + 42, x, domain='ZZ') Poly(x + 79, x, > domain='ZZ')⎥ > ⎢ > ⎥ > ⎢Poly(x + 98, x, domain='ZZ') Poly(x + 84, x, domain='ZZ') Poly(x + > 54, x, domain='ZZ') Poly(x + 32, x, domain='ZZ') Poly(x + 52, x, > domain='ZZ')⎥ > ⎢ > ⎥ > ⎣Poly(x + 10, x, domain='ZZ') Poly(x + 87, x, domain='ZZ') Poly(x + > 27, x, domain='ZZ') Poly(x + 23, x, domain='ZZ') Poly(x + 52, x, > domain='ZZ')⎦ > > In [40]: %timeit D = C.T * C > 100 loops, best of 3: 5.72 ms per loop > > On May 29, 3:44 pm, SherjilOzair <sherjiloz...@gmail.com> wrote: > > > > > > > > > Thanks for the pointers. I'm working on what you said. > > > Here are a few questions. > > > 1. Poly(123) should not give an error, why not treat it as a constant > > polynomial ? > > > 2. I used construct_domain on the list of elements of the Matrix, It > > returned DMP type, which does *not* allow coercion. What to do if one > > of my algorithms involve a square root ? > > > 3. Even when I tried operating on DMPs using an algorithm which did > > not have a square root, an "unsupported operand type(s) for /: 'int' > > and 'DMP'" TypeError was returned. > > > 4. Should I write different code for the algorithms for each > > groundtype ? For example, when using Sympy's type, using Add(*(...)) > > adds efficiency, but it doesn't make sense for other types. I can use > > sum(...) but that will sacrifice performance slightly. > > > 5. I think coercion is important for matrix algorithms. Other than > > Poly, which other lower-level classes allow coercion ? > > > - Sherjil Ozair > > > On May 28, 11:34 pm, Mateusz Paprocki <matt...@gmail.com> wrote: > > > > Hi, > > > > On 28 May 2011 15:30, SherjilOzair <sherjiloz...@gmail.com> 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#benc... > > > > > 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 sympy@googlegroups.com. > > > > To unsubscribe from this group, send email to > > > > sympy+unsubscr...@googlegroups.com. > > > > 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 sympy@googlegroups.com. To unsubscribe from this group, send email to sympy+unsubscr...@googlegroups.com. For more options, visit this group at http://groups.google.com/group/sympy?hl=en.