Speeding up multiplication in a free algebra.

[broken_symlink]

I've done some profiling of my code and it seems that my main slowdowns are coming from multiplication and addition of free algebra elements. Is there anyway I can make it go faster? I would like to do some more complicated computations but currently my code has been running for half an hour and still hasn't produced any results. I plan on parallelizing it in a few places, but was wondering if I could do anything about this as well. I tried compiling phi_ext_sk down to c with cython and added type hints in a few places but it didn't help. 

Specifically, here is what I get.

   307                                                   """

   308                                                   phi_b_ext returns phi_b_ext(a_ij)

   309                                                   """

   310                                           

   311        32        37956   1186.1      0.2          e = self.phi_ext_sk(self.w[0], a)

   312                                           

   313       416          431      1.0      0.0          for i in self.w[1:]:

   314       384          229      0.6      0.0              new_e = 0

   315     10152        11593      1.1      0.1              for k, v in e._FreeAlgebraElement__monomial_coefficients.iteritems():

   316      9768        14406      1.5      0.1                  new_g = 1

   317     58800       158459      2.7      0.7                  for g, p in k:

   318     49032      6465240    131.9     30.0                      g1 = self.phi_ext_sk(i, g)

   319     49032      7820633    159.5     36.3                      new_g = new_g*g1

   320     49032        47335      1.0      0.2                      if p != 1:

   321                                                                                 print "POW!!!!"

   322      9768      2818386    288.5     13.1                  new_g = v*new_g

   323      9768      4146233    424.5     19.3                  new_e = new_e + new_g

   324       384         6524     17.0      0.0              e = new_e

   325                                           

   326        32           16      0.5      0.0          return e

Function: phi_ext_sk at line 253                                                                                                                                       [4/681]

Total time: 6.70627 s

Line #      Hits         Time  Per Hit   % Time  Line Contents

==============================================================

   253                                               @profile

   254                                               def phi_ext_sk(self, k, a):

   255                                                   """

   256                                                   phi_ext_sk returns phi^ext_sigma_z(a_xy) 

   257                                                   """

   258                                           

   259     49064      3622876     73.8     54.0          i, j = self.gdicta[a]

   260                                           

   261     49064        46425      0.9      0.7          if k > 0:

   262                                           

   263     49064        34162      0.7      0.5              if i == k + 1:

   264     12004         8191      0.7      0.1                  if j != k and j != k + 1:

   265      8460         8997      1.1      0.1                      return self.astar[k, j]

   266                                                           else:

   267      3544         3904      1.1      0.1                      return self.astar[k, k+1]

   268     37060        23978      0.6      0.4              elif j == k + 1:

   269     12004         8038      0.7      0.1                  if i != k and i != k + 1:

   270      8460         9035      1.1      0.1                      return self.astar[i, k]

   271                                                           else:

   272      3544         3927      1.1      0.1                      return self.astar[k+1, k]

   273     25056        16273      0.6      0.2              elif i == k and j != k and j != k + 1:

   274      7016       685746     97.7     10.2                  return -1*self.astar[k+1, j]-self.astar[k+1, k]* \

   275      7016       760791    108.4     11.3                      self.astar[k, j]

   276     18040        12076      0.7      0.2              elif j == k and i != k and i != k + 1:

   277      7016       681617     97.2     10.2                  return -1*self.astar[i, k+1]-self.astar[i, k]* \

   278      7016       759046    108.2     11.3                      self.astar[k, k+1]

   279     11024         9371      0.9      0.1              elif i != k and i != k+1 and j != k and j != k+1:

   280     11024        11820      1.1      0.2                  return self.astar[i, j]

https://github.com/syamajala/KnotContactHomology/blob/master/braids.py

[Simon King]

Hi,

On 2013-07-10, broken_symlink <syam...@gmail.com> wrote:
> I've done some profiling of my code and it seems that my main slowdowns are
> coming from multiplication and addition of free algebra elements.
> Is there
> anyway I can make it go faster?

You define

self.alg = FreeAlgebra(self.br, len(g.split(',')), g)

This yields a very generic and very slow implementation, in particular
by the method of constructing elements later in your code (see below).

If all your elements are homogeneous (or at least: weighted homogeneous with respect
to positive integer weights of the algebra generators), then you might instead use

self.alg = FreeAlgebra(self.br, len(g.split(',')), g, implementation='letterplace')

I think their arithmetic is faster, and in particular you can do weighted
homogeneous Gröbner basis computations (which might be interesting in
your application, I don't know).

If your algebra elements are not weighted homogeneous, then you might
consider to use gap (via pexpect) or libgap (i.e., would be possible to
do on C library level) as a back-end. Gap has an implementation of free
associative unital algebras, and I think it is faster than the native
implementation in Sage. However, the (lib)gap free algebras are not
wrapped in Sage yet, hence, you would need to speak with (lib)gap on a
lower level.

One remark:

In your code, you construct an algebra element by
FreeAlgebraElement(self.alg, k)
Please don't do this!
In general, elements of an algebraic structure "A" given by data "d" should be
created via "A(d)" or via arithmetic operations on A.gens(), but by calling some
class constructor.

The reason is that A knows which implementation to use for the elements.
For example, if A happens to be a free algebra with Letterplace as back-end,
then the elements of A will *not* be instances of FreeAlgebraElement!
And if you use the generic (slow) implementation, then the elements will
use a *sub*class of FreeAlgebraElement. Example:

sage: FL = FreeAlgebra(QQ, 'a,b,c', 3, implementation='letterplace')
sage: type(FL.0).mro()
[sage.algebras.letterplace.free_algebra_element_letterplace.FreeAlgebraElement_letterplace,
sage.structure.element.AlgebraElement,
sage.structure.element.RingElement,
sage.structure.element.ModuleElement,
sage.structure.element.Element,
sage.structure.sage_object.SageObject,
object]
sage: FG = FreeAlgebra(QQ, 'a,b,c', 3)
sage: type(FG.0).mro()
[sage.algebras.free_algebra_element.FreeAlgebra_generic_with_category.element_class,
sage.algebras.free_algebra_element.FreeAlgebraElement,
sage.structure.element.AlgebraElement,
...
]

Best regards,
Simon

[Simon King]

PS:

On 2013-07-10, Simon King <simon...@uni-jena.de> wrote:
> One remark:
>
> In your code, you construct an algebra element by
> FreeAlgebraElement(self.alg, k)
> Please don't do this!
> In general, elements of an algebraic structure "A" given by data "d" should be
> created via "A(d)" or via arithmetic operations on A.gens(), but by calling some
> class constructor.

Oops, a typo: ..., but *not* by calling some class constructor.

Cheers,
Simon

[broken_symlink]

On Wednesday, July 10, 2013 5:20:29 AM UTC-4, Simon King wrote:

Hi,

On 2013-07-10, broken_symlink <syam...@gmail.com> wrote:
> I've done some profiling of my code and it seems that my main slowdowns are
> coming from multiplication and addition of free algebra elements.
> Is there
> anyway I can make it go faster?

You define

  self.alg = FreeAlgebra(self.br, len(g.split(',')), g)

This yields a very generic and very slow implementation, in particular
by the method of constructing elements later in your code (see below).

If all your elements are homogeneous (or at least: weighted homogeneous with respect
to positive integer weights of the algebra generators), then you might instead use

  self.alg = FreeAlgebra(self.br, len(g.split(',')), g, implementation='letterplace')

I think their arithmetic is faster, and in particular you can do weighted
homogeneous Gröbner basis computations (which might be interesting in
your application, I don't know).

My elements are not homogeneous.

 

If your algebra elements are not weighted homogeneous, then you might
consider to use gap (via pexpect) or libgap (i.e., would be possible to
do on C library level) as a back-end. Gap has an implementation of free
associative unital algebras, and I think it is faster than the native
implementation in Sage. However, the (lib)gap free algebras are not
wrapped in Sage yet, hence, you would need to speak with (lib)gap on a
lower level. 

I'll look into this. It seems like it would require a large rewrite though.

 

 

One remark:

In your code, you construct an algebra element by
  FreeAlgebraElement(self.alg, k)
Please don't do this!
In general, elements of an algebraic structure "A" given by data "d" should be
created via "A(d)" or via arithmetic operations on A.gens(), but by calling some
class constructor.

Thanks! I fixed that by doing for g, p in k._element_list and calling self.alg.gen(g)**p also i didn't realize I wasn't taking the exponent into account.  

On another note, I was able to get good speedups by parallelizing some parts, I don't know how long it will take to do the computation I'm interested in though as I haven't tried yet. 

[broken_symlink]

Is the code it uses for multiplication in free_algebra_element.py under devel/sage/sage/algebras? If so, I'm pretty sure for multiplication at least it should be possible to do better. Is there a way I can mess around with the _mul_ function easily? 

[Simon King]

Hi!

On 2013-07-14, broken_symlink <syam...@gmail.com> wrote:
> Is the code it uses for multiplication in free_algebra_element.py under
> devel/sage/sage/algebras? If so, I'm pretty sure for multiplication at
> least it should be possible to do better. Is there a way I can mess around
> with the _mul_ function easily?

I think so. If you want to be sure: Create an element x of a free algebra,
and then do
sage: edit(x._mul_, 'vim')
(where 'vim' can be replaced by the name of your favourite editor), and
then you can simply edit the sources.

Of course, you could also use inspection to find out the source file and
edit it outside of a sage session:

sage: F.<a,b> = FreeAlgebra(QQ)
sage: import inspect
sage: inspect.getfile(a._mul_)
'.../local/lib/python2.7/site-packages/sage/algebras/free_algebra_element.py'

Note: Do *not* edit the __mul__ (double underscore) method,
which is code that is (or should be) common for all elements in
Sage.
And also note: For Cython code, the "inspect" module would not always
give the right answer. But you can use
sage.misc.sageinspect.sage_getfile and similar methods, and of course
you could simply inspect the source code by

sage: a._mul_??

Anyway, when you have edited the _mul_ (*single* underscore) method, you
can test it after doing "sage -b" or "sage -br" (the latter will start a
new Sage session after rebuilding the changes).

And if you are happy with your changes, you can get a trac account, open
a ticket, post a patch there, and get someone to review your changes.
More details can be found in the "Sage developers guide".

Best regards,
Simon

[broken_symlink]

Thanks!

I have something that is passing the tests. I thought I would be able to do better by getting rid of the two for loops and using list comprehensions. But, according to kernprof all I managed to do was make things worse. 

    def _mul_(self, y):                                                               

        """                                                                           

        Return product of self and y (another free algebra element with the           

        same parents)                                                                 

                                                                                      

        EXAMPLES::                                                                    

                                                                                      

            sage: A.<x,y,z>=FreeAlgebra(ZZ,3)                                         

            sage: (x+y+x*y)*(x+y+1)    # indirect doctest                             

            x + y + x^2 + 2*x*y + y*x + y^2 + x*y*x + x*y^2                           

        """                                                                           

                                                                                      

        A = self.parent()                                                             

        z_elt = {}                                                                    

                                                                                      

        keys = product(self.__monomial_coefficients.keys(), y.__monomial_coefficients\

.keys())                                                                              

        keys = [reduce(lambda x, y: x*y, i) for i in keys]                            

                                                                                      

        values = product(self.__monomial_coefficients.values(), y.__monomial_coeffici\

ents.values())                                                                        

        values = [reduce(lambda x, y: x*y, i) for i in values]                        

                                                                                      

        for i in range(len(keys)):                                                    

            key = keys[i]                                                             

            if key in z_elt:                                                          

                z_elt[key] += values[i]                                               

            else:                                                                     

                z_elt[key] = values[i]  

        z = A(0)                                                                      

        z.__monomial_coefficients = z_elt                                             

        return z   

[Simon King]

Hi!

On 2013-07-14, broken_symlink <syam...@gmail.com> wrote:
> I have something that is passing the tests. I thought I would be able to do
> better by getting rid of the two for loops and using list comprehensions.
> But, according to kernprof all I managed to do was make things worse.

Anyway, if at some point you obtain an improvement: The right location
to post code is the Sage trac.

Best regards,
Simon

[broken_symlink]

I tried a lot of different things but was not able to do any better. I ended up just implementing my program using a different algorithm. It doesnt seem to be working in all cases yet, but its 2x faster than the old program in the cases its producing the same answer as the old program, and i haven't even done any optimizations yet. Thanks for all the help.