Opened 3 years ago

Closed 3 years ago

#26080 closed enhancement (fixed)

The Baker-Campbell-Hausdorff formula for nilpotent Lie algebras

Reported by: gh-ehaka Owned by:
Priority: major Milestone: sage-8.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:

Status badges

Description

Computation of the Baker-Campbell-Hausdorff 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 gh-ehaka

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.

comment:2 Changed 3 years ago by gh-ehaka

  • Keywords Lie algebras BCH formula added

comment:3 Changed 3 years ago by gh-ehaka

  • Authors set to Eero Hakavuori
  • Branch set to u/gh-ehaka/bch-26080
  • 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 Casas-Murua (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:

1276e0etrac #26080: initial implementation of the BCH formula

comment:4 Changed 3 years ago by tscrim

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 stress-tested 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., using yield instead of return).
  • 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/last-names (e.g., it is pbw_basis, not PBW_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 git

  • Commit changed from 1276e0eb34fc6509ef523cedee68803968fa70aa to 5972b6b0fbe210019e72e844a1420c7455a7f03f

Branch pushed to git repo; I updated commit sha1. New commits:

5972b6btrac #26080: converted bch iterator to a generator function and added interface for non-nilpotent Lie algebras

comment:6 Changed 3 years ago by gh-ehaka

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 tscrim

  • Reviewers set to Travis Scrimshaw

We are all always learning. :)

Four last little things I missed on the first go-around (but this is it, promise):

  • (m - 1) / 2 + 1 -> (m - 1) // 2 + 1 because on Python3, that will become a floating point and range 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 for coeff.
  • You can use the lower level IntegerListsLex to iterate over:
    sage: IV = IntegerVectorsConstraints(10 - 1, 2 * 3, min_part=1)
    sage: ILL = IntegerListsLex(10-1, 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 of Element, which have extra overhead from not being a builtin Python type.
  • It is faster to use X.bracket(Y) instead of L.bracket(X, Y) as the latter calls the former. If you actually know the elements belong to the same parent, you can do X._bracket_(Y) and avoid some extra overhead checks, but that should be negligible.

comment:8 Changed 3 years ago by git

  • Commit changed from 5972b6b0fbe210019e72e844a1420c7455a7f03f to e12ea13e8fc20f818b6b13afd604e5fb8913f793

Branch pushed to git repo; I updated commit sha1. New commits:

e12ea13trac #26080: efficiency improvements

comment:9 Changed 3 years ago by gh-ehaka

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 tscrim

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.

Last edited 3 years ago by tscrim (previous) (diff)

comment:11 follow-up: Changed 3 years ago by gh-ehaka

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 ; follow-up: Changed 3 years ago by tscrim

Replying to gh-ehaka:

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 finite-dimensional 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 gh-ehaka

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 tscrim

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 git

  • Commit changed from e12ea13e8fc20f818b6b13afd604e5fb8913f793 to 369f7a626c6aa4e362ac3d671ca4e28963d9163d

Branch pushed to git repo; I updated commit sha1. New commits:

369f7a6trac #26080: replaced old stopiteration

comment:16 Changed 3 years ago by gh-ehaka

At least with the current recursion algorithm, it is possible that all of the brackets [X-Y,Z_{m-1}] and [Z_{k_1},...[Z_{k_{2p}}, X+Y]...] computed in the recursion are zero for some m, but nonetheless higher order terms are non-zero. 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 [X-Y,Z_2], which has been quotiented out. Yet there are non-zero 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 tscrim

  • 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 vbraun

  • Branch changed from u/gh-ehaka/bch-26080 to 369f7a626c6aa4e362ac3d671ca4e28963d9163d
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.