Oh. I understand. And how does one turn off cahe ? Sent on my BlackBerry® from Vodafone
-----Original Message----- From: Mateusz Paprocki <matt...@gmail.com> Sender: sympy@googlegroups.com Date: Sun, 29 May 2011 22:37:30 To: <sympy@googlegroups.com> Reply-To: sympy@googlegroups.com Subject: Re: [sympy] Re: Borrowing ideas from Polys to Matrix Hi, On 29 May 2011 12:49, 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 some cases overhead of Poly maybe significant. But I'm quite sure you just didn't turn off cache, because for example: In [4]: f = x In [5]: g = Poly(x) In [6]: %timeit u = f + f 10000 loops, best of 3: 114 us per loop In [7]: %timeit u = g + g 10000 loops, best of 3: 40.2 us per loop Your example gives now: In [8]: A = randMatrix(5,5) In [9]: A = A.applyfunc(lambda i: i + x) In [10]: %timeit B = A.T * A 10 loops, best of 3: 20.6 ms per loop In [11]: C = A.applyfunc(poly) In [12]: %timeit D = C.T * C 100 loops, best of 3: 11.9 ms per loop So Poly is 2x faster than Add in this case. If you use %timeit and don't turn off cache then what you time is how fast @cache can retrieve previous results. > > 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. > > 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. -- 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.
