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.
