Oh. I understand. And how does one turn off cahe ?
Sent on my BlackBerry® from Vodafone
------------------------------------------------------------------------
*From: * Mateusz Paprocki <[email protected]>
*Sender: * [email protected]
*Date: *Sun, 29 May 2011 22:37:30 +0200
*To: *<[email protected]>
*ReplyTo: * [email protected]
*Subject: *Re: [sympy] Re: Borrowing ideas from Polys to Matrix
Hi,
On 29 May 2011 12:49, SherjilOzair <[email protected]
<mailto:[email protected]>> 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 <[email protected]
<mailto:[email protected]>> 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 <[email protected]
<mailto:[email protected]>> wrote:
>
>
>
>
>
>
>
> > Hi,
>
> > On 28 May 2011 15:30, SherjilOzair <[email protected]
<mailto:[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#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 [email protected]
<mailto:[email protected]>.
> > > To unsubscribe from this group, send email to
> > > [email protected]
<mailto:sympy%[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]
<mailto:[email protected]>.
To unsubscribe from this group, send email to
[email protected]
<mailto:sympy%[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.
--
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.
