Opened 3 years ago
Closed 3 years ago
#26080 closed enhancement (fixed)
The BakerCampbellHausdorff formula for nilpotent Lie algebras
Reported by:  ghehaka  Owned by:  

Priority:  major  Milestone:  sage8.4 
Component:  algebra  Keywords:  Lie algebras, BCH formula 
Cc:  Merged in:  
Authors:  Eero Hakavuori  Reviewers:  Travis Scrimshaw 
Report Upstream:  N/A  Work issues:  
Branch:  369f7a6 (Commits, GitHub, GitLab)  Commit:  369f7a626c6aa4e362ac3d671ca4e28963d9163d 
Dependencies:  #26076  Stopgaps: 
Description
Computation of the BakerCampbellHausdorff formula Z = log(exp(X)exp(Y))
in the case of nilpotent Lie algebras, where it reduces to a (complicated) polynomial expression.
Change History (18)
comment:1 Changed 3 years ago by
comment:2 Changed 3 years ago by
 Keywords Lie algebras BCH formula added
comment:3 Changed 3 years ago by
 Branch set to u/ghehaka/bch26080
 Commit set to 1276e0eb34fc6509ef523cedee68803968fa70aa
 Dependencies changed from #26076, #26036 to #26076
 Status changed from new to needs_review
An initial implementation is in the commits. I ended up making the computation into an iterator returning the terms involving brackets of increasing lengths, so it is in theory possible to compute the BCH formula to arbitrary precision in any Lie algebra. I added a BCH
method to LieAlgebras.Nilpotent
to compute the full expression using this iterator.
The iterator might be of independent interest for the computation of the (low order terms of the) abstract BCH formula in a free Lie algebra, however I did not know what would be the best way to expose it to the user.
One unfortunate part is that the computation time becomes quite unwieldy for the later terms. For example on my computer, computing the terms of increasing degree for the free Lie algebra on 2 generators took the following times:
10 : 45.4 ms 11 : 116 ms 12 : 327 ms 13 : 1.08 s 14 : 4.2 s 15 : 16.3 s 16 : 1min 8s 17 : 4min 59s 18 : 20min 52s
The only more efficient algorithm I know of is the one of CasasMurua (https://arxiv.org/abs/0810.2656). Based on their statement of requiring <15 mins for all terms of degree up to 20, I would guess that the efficiency improvement would be somewhere in the order of a 40x speed increase. For now, I think the simpler but less efficient approach is sufficient though.
New commits:
1276e0e  trac #26080: initial implementation of the BCH formula

comment:4 Changed 3 years ago by
Probably a good part of it might be a poor implementation of the Lie bracket for the free Lie algebra. At least it is taking on one iteration 14s/19s. Some caching, better implementation, and more Cythonization might help. At least, how it is currently designed is not very smart. For instance, it might be better to store the basis elements actually as lists that keeps track of the tree structure through some other list. There are definitely ways it can be improved (this is the first time that code is really being stresstested for efficiency).
Anyways, I think this will work for now as we can always work on it more on a followup. Some comments on the code.
 I don't see the reason why the
bch_formula
should be a class. I think it would be better suited as just a generator function (i.e., usingyield
instead ofreturn
).  Break the really long
.. MATH::
line. Also avoid using\dots
and instead be more specific about what kind of dots (in this case,\cdots
) (I don't really trust\dots
).  You should have an alias
baker_campbell_hausdorff = BCH
. Also, Sage convention is still method names are generally lower cases, even for acronyms/lastnames (e.g., it ispbw_basis
, notPBW_basis
).  Instead of just having a method for nilpotent Lie algebras, I think it would be good to have one that allows for a "precision" to be inputted. For instance
def bch(self, X, Y, prec=None): if self not in LieAlgebras.Nilpotent and prec is None: raise ValueError("the Lie algebra is not known to be nilpotent," " so you must specify the precision") from sage.algebras.lie_algebras.bch import BCH_iterator if prec is None: return self.sum(Z for Z in BCH_iterator(X, Y)) bch = BCH_iterator(X, Y) return self.sum(next(bch) for k in range(prec)) # maybe prec+1?
comment:5 Changed 3 years ago by
 Commit changed from 1276e0eb34fc6509ef523cedee68803968fa70aa to 5972b6b0fbe210019e72e844a1420c7455a7f03f
Branch pushed to git repo; I updated commit sha1. New commits:
5972b6b  trac #26080: converted bch iterator to a generator function and added interface for nonnilpotent Lie algebras

comment:6 Changed 3 years ago by
Thanks for the suggestions, changing the iterator to a generator seemed to indeed make a lot of sense since it cut out a lot of the unnecessary pieces of code! Also good to know about the naming convention, I wasn't sure what the policy was for names.
comment:7 Changed 3 years ago by
 Reviewers set to Travis Scrimshaw
We are all always learning. :)
Four last little things I missed on the first goaround (but this is it, promise):
(m  1) / 2 + 1
>(m  1) // 2 + 1
because on Python3, that will become a floating point andrange
does not like floating point limits (//
is floor division). A faster way to create fractions in
QQ
:sage: %timeit QQ(1) / QQ(121) The slowest run took 12.50 times longer than the fastest. This could mean that an intermediate result is being cached. 1000000 loops, best of 3: 1.11 µs per loop sage: %timeit QQ((1, 121)) The slowest run took 21.95 times longer than the fastest. This could mean that an intermediate result is being cached. 1000000 loops, best of 3: 771 ns per loop sage: %timeit ~QQ(121) The slowest run took 27.61 times longer than the fastest. This could mean that an intermediate result is being cached. 1000000 loops, best of 3: 440 ns per loop
Similarly forcoeff
.  You can use the lower level
IntegerListsLex
to iterate over:sage: IV = IntegerVectorsConstraints(10  1, 2 * 3, min_part=1) sage: ILL = IntegerListsLex(101, length=2*3, min_part=1) sage: [list(x) for x in IV] == [list(x) for x in ILL] True
IntegerVectorsConstraints
is slower because it generates instances ofElement
, which have extra overhead from not being a builtin Python type.  It is faster to use
X.bracket(Y)
instead ofL.bracket(X, Y)
as the latter calls the former. If you actually know the elements belong to the same parent, you can doX._bracket_(Y)
and avoid some extra overhead checks, but that should be negligible.
comment:8 Changed 3 years ago by
 Commit changed from 5972b6b0fbe210019e72e844a1420c7455a7f03f to e12ea13e8fc20f818b6b13afd604e5fb8913f793
Branch pushed to git repo; I updated commit sha1. New commits:
e12ea13  trac #26080: efficiency improvements

comment:9 Changed 3 years ago by
I only now realized that bernoulli
already returns elements of QQ
so the first cast in coeff
was unnecessary, so I left the latter cast and division in place. The other improvements are now in.
comment:10 Changed 3 years ago by
Thanks. I unintentionally lied: raise StopIteration
> return
. You can set a positive review on my behalf once fixed.
Something we might want to consider for a followup ticket (you can do it here if you want) is when, say, a Lie algebra does not know it is nilpotent or X,Y
generate a finite BCH formula, we could have a test to stop the iteration. I think this would be just checking that all of the brackets computed at a particular depth are 0. I leave it up to you.
comment:11 followup: ↓ 12 Changed 3 years ago by
That sounds like a good idea to add. It would also make sense then to allow prec=Infinity
to handle the usecase where the formula is known to be finite even if the algebra is not nilpotent. Is there a canonical format for an input that can be an integer or infinite?
I'll try to edit it into this ticket later today or tomorrow.
comment:12 in reply to: ↑ 11 ; followup: ↓ 13 Changed 3 years ago by
Replying to ghehaka:
That sounds like a good idea to add. It would also make sense then to allow
prec=Infinity
to handle the usecase where the formula is known to be finite even if the algebra is not nilpotent.
The problem with that there is no way to test whether a Lie algebra is nilpotent or not in general (at least, as far as I know). For the finitedimensional ones, we can compute the lower central series of course, but that is expensive. So IMO it is better to avoid the computation and to play it safe. However, that is a very weak position as I can see good reason to just do the check in that case. Your choice.
Is there a canonical format for an input that can be an integer or infinite?
Not to my knowledge. You have to check both x in ZZ or x == Infinity
.
comment:13 in reply to: ↑ 12 Changed 3 years ago by
I meant the usecase when the user knows for some a priori reason that the formula is finite, but maybe does not know the maximum degree. So there would be no nilpotency check, but the program would just compute until the formula terminates.
This would of course mean that if passed 'bad' parameters, the method would hang. Would this be bad practice to allow even with a warning about it in the doc?
Although thinking about it a bit more, it is no longer clear to me if all brackets being zero on a degree is sufficient for termination. The brackets in the recursion are not very transparent, so it is not clear to me if there could be some corner case that fails. I will have to think about it more, and if I can't figure it out in a reasonable time, it can be left to a later improvement :)
comment:14 Changed 3 years ago by
Ah, I see what you mean. Yes, having prec=oo
to let it run is fine. As you said, we just put a warning in the doc that it may run forever, but the user then knows what s/he is doing. +1 for this.
comment:15 Changed 3 years ago by
 Commit changed from e12ea13e8fc20f818b6b13afd604e5fb8913f793 to 369f7a626c6aa4e362ac3d671ca4e28963d9163d
Branch pushed to git repo; I updated commit sha1. New commits:
369f7a6  trac #26080: replaced old stopiteration

comment:16 Changed 3 years ago by
At least with the current recursion algorithm, it is possible that all of the brackets [XY,Z_{m1}]
and [Z_{k_1},...[Z_{k_{2p}}, X+Y]...]
computed in the recursion are zero for some m
, but nonetheless higher order terms are nonzero. An example of this is given by quotienting a Lie algebra by the term of degree 3 of the BCH formula:
sage: L = LieAlgebra(QQ, 2, step=4) sage: X, Y = L.homogeneous_component_basis(3) sage: Q = L.quotient(X + Y) sage: Q.inject_variables() Defining X_1, X_2, X_12, X_112, X_1112 sage: Q.bch(X_1, X_2) X_1 + X_2 + 1/2*X_12  1/24*X_1112
For third order terms, the IntegerListsLex
part of the sum is only [Z_1, [Z_1, Z_1]] = 0
, so Z_3
will be just a multiple of [XY,Z_2]
, which has been quotiented out. Yet there are nonzero terms of degree 4.
At this point I don't know what would be a possible alternative check to see if the formula terminates (other than testing for nilpotency of the Lie algebra generated by the two elements). So I would say to leave the improvements to a later ticket.
comment:17 Changed 3 years ago by
 Status changed from needs_review to positive_review
This will definitely work for now (and IIRC, it sufficient for your purposes). Thank you for looking into it (and that interesting example).
comment:18 Changed 3 years ago by
 Branch changed from u/ghehaka/bch26080 to 369f7a626c6aa4e362ac3d671ca4e28963d9163d
 Resolution set to fixed
 Status changed from positive_review to closed
Some implementation already exists in a codebase that I need to clean up and import into Sage. The method of computing is to compute the BCH formula for 2 generators of a free nilpotent Lie algebra of the same nilpotency step and to map the result to the desired nilpotent Lie algebra via a Lie algebra morphism mapping the generators to X and Y.