The iisignature library: efficient calculation of iterated-integral signatures and log signatures

Introduction

iisignature is a Python package which calculates the iterated-integral signatures and log signatures of paths. The signature is an object which is crucial in the mathematical theory of rough paths, and the calculations have proved to be useful in machine-learning applications, particularly classification problems where the data itself is a stream or a path in space, ranging from an application to online Chinese handwriting recognition in 2013 [BEN] to skeleton-based human action recognition in 2017 [action]. Other domains where the data has this form include signals from EEG and other medical monitors, sound and financial time series, where some set of numbers is varying in time. Often the samples can be noisy, can have varying length and both local and global structure can be important. A survey of such applications is given in [OxSigIntro].

An existing open-source implementation is the esig package from CoRoPa[coropa]. CoRoPa operates in a sparse fashion, keeping track of only non-zero elements of the signature. This has been known as sparse signatures. It is useful in some applications of signatures in high-dimensional spaces where the path only moves in certain combinations of the input dimensions.

The particular focus of iisignature is piecewise-linear fixed-dimensional paths which typically move in all their dimensions. In this setting, usually none of the elements of the signature are zero. Sparse methods impose a significant overhead in this context; iisignature is directed at these dense signatures.

We study the mathematical properties of the free Lie algebra to implement algorithms for calculating signatures in the dense case. We also benchmark the performance of these algorithms, and provide an efficient open-source implementation. This paper is organised as follows. The rest of this section introduces signatures, log signatures and the library. Signature algorithms are discussed in 2. Log signature methods are introduced in 3, the direct method is discussed in 4 and the projection method in 5. Considerations around the implementation are presented in . Indicative timings are given in and memory usages in . We briefly discuss other functionality provided by the library in before concluding.

iisignature is hosted at https://github.com/bottler/iisignature and is available on PyPI.

The iterated-integral signature of a continuous path is an infinite sequence of numbers. It is used in the mathematical theory of differential equations driven by paths. In these problems, a path is the driving signal for a certain type of system. It turns out that the signature is exactly the information about a path which you need to know in order to predict how the output of the system will behave, using a generalisation of Taylor’s theorem. It is natural that the signature would also be the right information to extract from a path if we want a machine-learning algorithm to understand the shape of the path.

The signature is divided into units called levels. We cannot store the whole signature of a path on a computer, rather we calculate a certain number of levels of it. The more levels of a signature are known, the more precisely the shape of the path is determined. If a path changes very slightly, the first few levels of its signature will also only change very slightly. If a path is moved (translated) but retains its shape, its signature will not change.

The number of elements of level $m$ of the signature of a $d$-dimensional path is $d^{m}$. They are the values of iterated integrals which consist of $m$ nested integrals, and they are labelled with $m$ numbers each corresponding to one of the dimensions. To distinguish these numbers which label the dimensions from other numbers, we write them bold and in blue. For example, a two-dimensional path might be given in coordinates as (\gamma_{{\color[rgb]{0,0,1}\mathbf{1}}}(t),\gamma_{{\color[rgb]{0,0,1}\mathbf{2}}}(t)) as $t$ varies from $a$ to $b$. Its signature is a function denoted by $X^{\gamma}_{a,b}$. Level three of its signature has eight elements, called X_{a,b}^{\gamma}({\color[rgb]{0,0,1}\mathbf{111}}), X_{a,b}^{\gamma}({\color[rgb]{0,0,1}\mathbf{112}}) and so on. The one indexed by the word {\color[rgb]{0,0,1}\mathbf{122}} is

The information in the first level of the signature is the total displacement of the path, i.e. the direction and distance from its starting point to its ending point. The information which the second level of the signature adds is the signed area of the path projected in each plane. Higher levels of the signature provide more detailed information about the path’s shape.

2 Signed area

For a two-dimensional path, the information carried by the first two levels of the signature is the total displacement of the path (in the first level, which is two numbers) and the signed area between the path and the straight line from its beginning to end. Figure 1 shows this information for two straight lines and their combination, which contains area.

The following is an intuitive definition of the signed area of a path in the plane. For a closed path, that is one which ends where it starts, the signed area is the sum of the signed areas of the regions bounded by the path, which is the area times the number of times the path goes round that region in an anticlockwise manner minus the number of times the path goes round it clockwise (i.e. the winding number). For example, in the path shown in Figure 2(a), regions whose areas count positively are labelled with a $+$, and negatively with a $-$. One region’s area counts twice negatively; it is labelled with $--$. For a more general path, its signed area is the signed area of the closed path you get by joining it with a straight line from its end to its start.

As an example of how the area can be useful in classifying the shape of the path, consider classifying handwritten digits 0 and 8. Usually these are written with a single stroke which ends near its beginning, so the displacement is insufficient for distinguishing them. However, the encompassed areas are statistically different. The figure 0 is typically formed from a single anticlockwise loop, generating a positive signed area, while the figure 8 contains two regions with opposite sign, leading to cancellation of signed area. The diagrams in Figure 2(b) and (c) illustrate this. The Pendigits dataset [pendigits] collected the traces of many people writing the digits 0 to 9, and the histogram in Figure 3 shows how different the signed areas of the first (and usually only) strokes of these digits are. This clear separation is an illustration of the potential usefulness of the signature for classification.

3 What is the log signature of a path?

The log signature is a compressed version of the signature. It carries the same information, but in a more compact way. It is also divided into levels. Up to level $m$, the log signature contains fewer numbers than the signature. Any given set of values for these numbers actually gives the log signature of some path, whereas this is not the case for signatures, because there is some redundancy in the signature. For example the first two levels of the signature of a two-dimensional path consists of $2+2^{2}=6$ numbers but we saw that this information is the path’s total displacement and signed area, which can be stored in three numbers, which are exactly the first two levels of the log signature. In applications, the log signature might be less susceptible to roundoff error. The log signature is defined in terms of the signature, in a way analogous to logarithms of numbers, but can be calculated via an independent algorithm.

4 Using the library

The library is designed to make calculating large numbers of signatures and log signatures fast. To this end, preparatory calculations for the log signature calculation happen in a separate preparation function called prepare. This also means the library’s size on disc can be small; there is no separate code for specific numbers of dimensions and levels.

Figure 4 shows an example of calculating the signature and log signature of a 3-dimensional path up to level 4, which is specified as a set of $n$ points.

After running this code, signature will be a numpy array of the values of levels 1, 2, 3 and 4 of the signature, which has length $3+3^{2}+3^{3}+3^{4}=120$. Note that level 0, which is the constant 1 and contains no information about the path, is excluded from the output of iisignature. logsignature will be a numpy array of the values of levels 1, 2, 3 and 4 of the logsignature, which has length $3+3+8+18=32$.

Signatures

Calculating the signature of a path can be done inductively relying on the following two rules.

If $\gamma$ is a straight line defined on the interval $[a,b]$ then its signature as a function on words is

Grouped by levels, using $x=\gamma(b)-\gamma(a)$ as the displacement, the signature looks like

where $\otimes$ is the tensor product. Alternatively, if each level is thought of as a vector of numbers, this formula should be read with $\otimes$ denoting the Kronecker product.

If $a<b<c$ then the result (from [chen]) known as Chen’s identity states that

Grouped by levels, this signature looks like

When calculating the signature of a path given as a series of straight-line displacements, we start with the signature of the first displacement (calculated from (2)) and step-by-step concatenate on the signature of each succeeding displacement using (4).

Level $m$ of the signature contains $d^{m}$ values. Calculating it for a displacement using (2) takes $d+d^{m}$ multiplications beyond what has already been calculated for lower levels. However, in the signature of a straight line, each level is a symmetric tensor and so level $m$ only contains $\binom{d+m-1}{m}$ distinct values, using the formula for unordered sampling with replacement. An alternative, more complicated, method that takes account of this redundancy exists. Only $d+\binom{d+m-1}{m}$ multiplications are required. Implementing it showed it to be slower, so iisignature does not use this idea.

Log Signatures

$\mathfrak{n}_{d,m}$ is a finite dimensional real vector space, but there is no single obvious basis for it. In order to use the log signature as an efficient representation of a path, we need to choose a fixed basis. There are two commonly used bases. They are both Hall bases[hall1950]. A Hall basis is made up of bracketed expressions, and it is determined by an ordering of all bracketed expressions.

The Lyndon basis[shirshov], which is the default in iisignature. Each basis element is labelled with a Lyndon word on \{{\color[rgb]{0,0,1}\mathbf{1}},{\color[rgb]{0,0,1}\mathbf{2}},\dots,{\color[rgb]{0,0,1}\mathbf{d}}\}, which is a sequence which comes earlier in lexicographic order than any of its rotations. (For example, the rotations of {\color[rgb]{0,0,1}\mathbf{2432}} are {\color[rgb]{0,0,1}\mathbf{2243}}, {\color[rgb]{0,0,1}\mathbf{3224}} and {\color[rgb]{0,0,1}\mathbf{4322}}. {\color[rgb]{0,0,1}\mathbf{2243}} and {\color[rgb]{0,0,1}\mathbf{1213}} are Lyndon words but {\color[rgb]{0,0,1}\mathbf{31}} and {\color[rgb]{0,0,1}\mathbf{3224}} are not.)

The standard/canonical Hall basis, which we implement in such a way as to match CoRoPa[coropa] exactly. The ordering of equal-length expressions $[A,B]$ and $[C,D]$ is defined recursively: $[A,B]<[C,D]$ if either $A<C$ or ($A=C$ and $B<D$).

In these bases, each basis element is either a letter or a single bracketed expression, whose left and right are basis elements. We always pick an order on basis elements such that shorter bracketed expressions come before longer ones, and single letters, which are the first level, are in their natural order {\color[rgb]{0,0,1}\mathbf{1}}<{\color[rgb]{0,0,1}\mathbf{2}}<\dots<{\color[rgb]{0,0,1}\mathbf{d}}.

Much of the algebra calculations can be done once in the prepare function. This is a major contribution of iisignature and ensures for given $d$, $m$ and the choice of basis that the calculation is as efficient as possible. This is relevant in machine learning applications where typically many similar calculations are required.

Log Signatures directly

The log signature of a straight line displacement is just the displacement itself in level 1, and zero in every other level. The log signature of the concatenation of two paths is the Baker-Campbell-Hausdorff (BCH) product of the log signatures of the two paths. The direct method for calculating the log signature relies on being able to transform the log signature of a path given in terms of one of the bases above to the log signature of that path concatenated with a fixed line segment, achieved using the BCH product.

The coefficients in this expansion up to terms of depth twenty have been calculated and distributed by Fernando Casas and Ander Murua at [bchinfo], using their method described in [bch]. We distribute their file as part of iisignature, and read it when necessary.

We can compute the Lie bracket of each pair of basis elements as a combination of other basis elements, and therefore, given two log signatures as combinations of basis elements (the second known to be just a displacement) we can find the expanded expression of their BCH product as a combination of basis elements. By doing this with indeterminates, the library develops an internal representation.

As an example, in the case where the Lyndon basis is used, and we are concerned with two dimensions up to level two, a log signature looks like

The inductive step of the algorithm to accumulate log signatures by adding linear segments for $d=m=2$ is shown in Figure 5.

If we go up to level 3, a log signature looks like

with the final algorithm being as shown in Figure 6.

These functions have a lot of common structure. First a sequence of monomials in the input elements are constructed in the temporary array $t$. Higher order monomials are calculated inductively from other elements of $t$ to deduplicate the necessary multiplications. Then some members of $a$ are incremented by some multiples of some of the temporary variables. Then the first $d$ elements of $a$ are incremented by all elements of $b$. Exactly which is given by the FunctionData structure. In general these functions are long and branching-free. The variable $a$ is modified in-place to produce the log signature of the extended path.

The basis (of the free Lie algebra on 2 symbols) used to express the BCH formula does not change the code we get, because the various equivalent bracketed expressions come to the same thing when they have been multiplied out. We use the Lyndon basis because it has slightly fewer terms, as [bch] describes and partially explains. This choice is independent of the choice of basis (of the free Lie algebra on $d$ symbols) in which the log signature is expressed. In general, we end up with fewer terms and a slightly faster calculation when the Lyndon basis is used for the log signature.

Log Signatures from Signatures

A simple method for calculating the log signature of a path is to calculate its signature first, and then convert to the log signature. The first step in doing the conversion is taking the logarithm itself in tensor space. This explicitly uses the formula (6) where $n$ only needs to go as high as the required level, and the power is in the concatenation product. This results in the log signature as an element of tensor space (which means it is as long as a signature), which is returned when logsig is called with the "X" (expanded) method. The exact order of evaluation of formula (6) for best efficiency which we use is one which was suggested by Mike Giles[Giles].

To express this Lie element into a specified basis, we need to project it. We calculate a projection explicitly. There are known explicit forms for projections, for example the map given by the Dynkin-Specht-Wever lemma directly ([DSWLemma]), which requires more operations. The prepare function calculates a projection upfront.

Given the bracketed expression of a basis element with $m$ letters, we can easily find its expression in expanded space, by multiplying out the brackets. For example, [[{\color[rgb]{0,0,1}\mathbf{1}},{\color[rgb]{0,0,1}\mathbf{3}}],{\color[rgb]{0,0,1}\mathbf{3}}] is {\color[rgb]{0,0,1}\mathbf{133}}-2\,{\color[rgb]{0,0,1}\mathbf{313}}+{\color[rgb]{0,0,1}\mathbf{331}}. This gives us the full matrix $M_{m}$ to transform each level of the log signature to its expanded version. Each column of $M_{m}$ is labelled with a basis element, and each row is labelled with one of the $d^{m}$ words of length $d$. To compress level $m$ a given expanded log signature $x_{m}$ to its value $c_{m}$ in terms of a basis, we just need to solve a least squares problem $M_{m}c_{m}=x_{m}$. This problem is a very overdetermined system which is known to have an exact answer, up to rounding considerations. $M_{m}$ is tall and skinny.

The words occurring in the terms of the expansion of such a bracketed expression are anagrams of the foliage of the expression. In that same example, for instance, {\color[rgb]{0,0,1}\mathbf{133}}, {\color[rgb]{0,0,1}\mathbf{313}} and {\color[rgb]{0,0,1}\mathbf{331}} are anagrams of {\color[rgb]{0,0,1}\mathbf{133}}. This leads to a lot of sparsity in the matrix $M_{m}$. Permuting the rows and columns to gather anagrams makes $M_{m}$ be a block diagonal matrix. We can save time doing the transformation by solving a separate linear system for each equivalence class of anagrams of words of length $m$.

For the standard Hall basis, this is exactly the procedure which we follow. In prepare, we determine all the mapping matrices between anagram classes of the log signature and its expansion, and then we calculate all their Moore-Penrose pseudoinverses, so that solving the systems is just a matrix multiplication. The number of words in an anagram set containing $m$ letters where the frequency of the $i$th letter is $n_{i}$ is given by a multinomial coefficient $\frac{m!}{n_{1}!\dots n_{d}!}$. The number of Lie basis elements in an anagram set is given by the second Witt formula of Satz 3 of [witt] as

where $\delta$ ranges over all common factors of the $n_{i}$ and $\mu$ is the Möbius function. In the simple special case that the words have $m$ distinct letters, there are $m!$ words and $(m-1)!$ basis elements. In the Lyndon case, this formula makes sense because the Lyndon words in such a set of $m!$ words are just all that begin with the lowest letter. Typically the largest anagram sets are the ones with about the same number of each letter. For them, (7) is just $\frac{1}{m}$ times the number of words in the set because 1 is the only value of $\delta$. For example, looking at level 10 for a 3-dimensional path, the signature has 59049 elements and the log signature 5880, and there are 63 anagram classes. The count is $63=\binom{10+3-1}{10}-3$ using the formula for unordered sampling with replacement and the fact that no basis element above level 1 has only one distinct letter in it. The 12 most balanced anagram classes account for 3708 elements of the log signature, or $63.1\%$ of it.

The big anagram classes account for most of the runtime when projecting to the log signature: multiplying a $420\times 4200$ matrix by a 4200-vector takes 80% more multiplications than multiplying a $312\times 3150$ matrix by a 3150-vector and so on.

If the Lyndon basis is required, then we have a more efficient implementation, which depends on a special property it has. On pages 89–91 of [FLA], the notation $P_{a}$ is introduced for the Lie polynomial corresponding to the Hall word $a$, i.e. the polynomial you get by multiplying out the bracketed expression corresponding to the unique basis element whose foliage is $a$. This notation is used in the statement of the following.

The simplest case of the final statement, where $w$ is itself a single Lyndon word, gives the following useful fact. When the bracketed expression corresponding to a Lyndon word is expanded and terms are collected and ordered in alphabetical order of the word, the first term will be the Lyndon word itself, with coefficient 1. (For an example, consider the Lyndon word {\color[rgb]{0,0,1}\mathbf{133}}; its bracketed expression is [[{\color[rgb]{0,0,1}\mathbf{1}},{\color[rgb]{0,0,1}\mathbf{3}}],{\color[rgb]{0,0,1}\mathbf{3}}] and we saw earlier that this expands to {\color[rgb]{0,0,1}\mathbf{133}}-2\,{\color[rgb]{0,0,1}\mathbf{313}}+{\color[rgb]{0,0,1}\mathbf{331}}.) This means that the tall skinny matrix $M_{m}$ is lower triangular, as are its anagram blocks. If we take such a block and remove all the rows corresponding to words which are not Lyndon, we are left with the mapping from an anagram class in the compressed log signature to same Lyndon word elements of the expanded signature. It is a square lower triangular matrix with ones on the diagonal. We can now solve the system directly in many fewer operations, with just addition and multiplication, just looking at the Lyndon word elements of the expanded signature. prepare determines the necessary indices and matrices, and logsig does the solving.

For example, in level 4 on 3 dimensions, the following are the three basis elements which contain two {\color[rgb]{0,0,1}\mathbf{1}}s, a {\color[rgb]{0,0,1}\mathbf{2}} and a {\color[rgb]{0,0,1}\mathbf{3}}:

The matrix corresponding to these looks as follows