---
title: 'pySigLib: Fast Signature Computation'
url: https://www.emergentmind.com/topics/pysiglib
type: topic
---

# pySigLib: Fast Signature Computation

Searching arXiv for the specified paper to ground the article.
pySigLib is a Python library for signature-based computation on sequential data, introduced as “Fast Signature-Based Computations on CPU and GPU” [2509.10613]. It provides optimised implementations of truncated signatures and signature kernels on CPU and GPU, is fully compatible with PyTorch’s automatic differentiation, and introduces a differentiation scheme for signature kernels that is described as exact, faster, and robust [2509.10613]. The library is positioned within the broader use of signature-based methods for time-series learning, especially where signature kernels are used as discriminators and training losses in generative modeling, including quantitative finance [2509.10613].

## 1. Conceptual setting and motivation

Signature-based methods map a path to a hierarchy of iterated integrals that linearize many nonlinear functionals on sequential data [2509.10613]. For a technical audience, the key point is that linear functionals of signatures are dense in continuous functionals on paths, up to reparameterization, which makes signatures a powerful representation for learning problems on streams and time series [2509.10613]. Signature kernels complement this viewpoint by providing inner products in the tensor algebra, thereby allowing access to the signature feature space implicitly and avoiding explicit construction of high-dimensional truncated feature maps [2509.10613].

The motivation for pySigLib is computational rather than conceptual. The paper states that prior libraries struggle to scale to the dataset sizes and sequence lengths encountered in practice [2509.10613]. The stated bottlenecks are threefold: direct signature computation grows rapidly with truncation degree and dimension; signature-kernel implementations often precompute refined grids and become memory-bound; and existing gradient computations for signature kernels rely on a second PDE whose discretization can yield inaccurate derivatives for short sequences or low refinement, harming training stability [2509.10613].

Within this setting, pySigLib is presented as an efficient software stack for large-scale signature-based computation. The paper attributes its performance to optimised CPU and GPU implementations, on-the-fly path transforms such as lead–lag and time augmentation, SIMD-friendly memory layouts, anti-diagonal GPU tiling, and batched matrix multiplication for kernel inner products [2509.10613]. A plausible implication is that the library is designed for use cases in which both throughput and gradient fidelity are central, particularly in end-to-end differentiable training pipelines.

## 2. Mathematical structure of signatures and signature kernels

For a \(d\)-dimensional path \(X : [0,T] \to \mathbb{R}^d\) of bounded variation, the truncated path signature up to degree \(M\) is defined as
$$
S(X) = (S^0(X), S^1(X), \ldots, S^M(X)),
$$
with
$$
S^m(X) = \int_{0 < t_1 < \cdots < t_m < T} dX_{t_1} \otimes \cdots \otimes dX_{t_m} \in (\mathbb{R}^d)^{\otimes m},
$$
and \(S^0(X)=1\) [2509.10613]. For piecewise linear interpolation of discrete samples \(X_{t_i}\), the paper identifies two algebraic facts as foundational for implementation: the signature of a linear segment with increment \(z = X_b - X_a\) is
$$
S(\text{segment}) = \sum_{k=0}^{\infty} \frac{z^{\otimes k}}{k!},
$$
and Chen’s identity for concatenation gives
$$
S(x * y)_{[a,c]} = S(x)_{[a,b]} \otimes S(y)_{[b,c]}.
$$
These identities underpin the library’s signature algorithms [2509.10613].

The signature kernel between two paths \(X\) and \(Y\) is the inner product in the tensor algebra. With weights \(\{\Lambda_m\}\), such as exponential weighting \(\Lambda_m = \lambda^m/m!\), the truncated kernel is
$$
k(X,Y) = \sum_{m=0}^{M} \Lambda_m \langle S^m(X), S^m(Y)\rangle,
$$
where the inner product is the natural one on tensor powers [2509.10613]. The paper emphasizes that pySigLib does not compute \(S(X)\) and \(S(Y)\) explicitly to obtain this kernel. Instead, it uses the Goursat PDE characterization
$$
\frac{\partial^2 K(s,t)}{\partial s \partial t} = \langle X'(s), Y'(t)\rangle K(s,t),
$$
with boundary conditions \(K(0,t)=K(s,0)=1\), and discretizes this PDE on dyadically refined grids [2509.10613].

For a practical discretization over grid indices \(i,j\), with increments
$$
\Delta_{i,j} = \langle X_{t_{i+1}} - X_{t_i}, Y_{s_{j+1}} - Y_{s_j}\rangle,
$$
the finite-difference update is
$$
\hat K_{i+1,j+1} = (\hat K_{i+1,j} + \hat K_{i,j+1}) A(\Delta_{i,j}) - \hat K_{i,j} B(\Delta_{i,j}),
$$
where
$$
A(\Delta) = 1 + \frac{1}{2}\Delta + \frac{1}{12}\Delta^2, \qquad
B(\Delta) = 1 - \frac{1}{12}\Delta^2.
$$
Dyadic refinement orders \((\lambda_1,\lambda_2)\) define a grid \(P_{\lambda_1,\lambda_2}\) of size \((2^{\lambda_1}L_X)\times(2^{\lambda_2}L_Y)\) [2509.10613]. This formulation is central to the kernel implementation and to the library’s backward pass.

## 3. Core algorithms and the differentiation scheme

pySigLib implements two principal algorithms for truncated signatures. The direct algorithm, described as iisignature-style, maintains the truncated signature \(A_0,\ldots,A_N\) in one flattened, contiguous array and updates each level in reverse order:
$$
A_k \leftarrow \sum_{i=0}^{k} A_i \otimes \frac{z^{\otimes (k-i)}}{(k-i)!}.
$$
Reverse-order updates permit in-place writes without temporary buffers, which minimizes allocations and memory traffic [2509.10613].

The second implementation is Horner’s method, described as signatory-style. The update is rearranged as
$$
A_k = (((z/k + A_1) \otimes z/(k-1) + A_2) \otimes z/(k-2) + \cdots \otimes z/2 + A_{k-1}) \otimes z + A_k.
$$
The implementation separates an intermediate
$$
B_k = (((z/k + A_1) \otimes z/(k-1) + A_2) \otimes \cdots \otimes z/2)
$$
and then performs \(B_k \otimes z + A_k\) in-place, reusing a single preallocated contiguous buffer sized for \(B_N\) and writing products in reverse order so that overwritten values are no longer needed [2509.10613]. The paper also states that backpropagation through signatures uses the time-reversed path to deconstruct the signature, applying analogous memory optimizations and a Horner-based deconstruction to reduce operations [2509.10613].

For signature kernels, the CPU implementation first precomputes
\[
dx = X_{1:L_X} - X_{0:L_X-1},
\]
and similarly \(dy\), then computes
\[
\Delta = dx^Tdy
\]
via batched matrix multiplication using `torch.bmm`; the paper states that this step dominates runtime for large \(d\) [2509.10613]. Dyadic refinement is then applied on-the-fly by referencing \(\Delta[s // 2^{\lambda_1}, t // 2^{\lambda_2}]\) at each stencil step, rather than precomputing refined grids [2509.10613].

The GPU implementation exploits the independence of entries on each anti-diagonal of the PDE grid. Only three anti-diagonals are kept and rotated in shared memory, and kernels are processed in blocks of 32 along the row or column dimension to avoid thread-count limits such as 1024 threads per block [2509.10613]. The first row is stored as an initial-condition vector of ones in global memory, and each block writes back the last row as the initial condition for the next block; with batches, blocks for different kernels run asynchronously [2509.10613].

The most distinctive algorithmic feature is the differentiation scheme for signature kernels. Existing methods are described as estimating gradients by solving a second PDE, which can be inaccurate for short sequences or low dyadic orders [2509.10613]. pySigLib instead differentiates the finite-difference solver itself to produce exact gradients in one backward pass over the same grid, reusing forward states [2509.10613]. For \(\lambda_1=\lambda_2=0\), the adjoint recurrences are given by
$$
\frac{\partial F}{\partial \hat K_{i,j}} =
\left(\frac{\partial F}{\partial \hat K_{i+1,j}}\right) A(\Delta_{i,j-1}) +
\left(\frac{\partial F}{\partial \hat K_{i,j+1}}\right) A(\Delta_{i-1,j}) -
\left(\frac{\partial F}{\partial \hat K_{i+1,j+1}}\right) B(\Delta_{i,j}),
$$
and
$$
\frac{\partial F}{\partial \Delta_{i,j}} =
\left(\frac{\partial F}{\partial \hat K_{i+1,j+1}}\right)
\left[
(\hat K_{i+1,j} + \hat K_{i,j+1})
\left(\frac{1}{2} + \frac{1}{6}\Delta_{i,j}\right)
+ \frac{1}{6}\hat K_{i,j}\Delta_{i,j}
\right].
$$
For dyadic refinement, multiple \(\hat K\) entries may depend on a given \(\Delta_{i,j}\), and the corresponding contributions are accumulated [2509.10613]. Gradients with respect to the increments satisfy
$$
\frac{\partial F}{\partial dx} = dy \left(\frac{\partial F}{\partial \Delta}\right)^T,\qquad
\frac{\partial F}{\partial dy} = dx \left(\frac{\partial F}{\partial \Delta}\right),
$$
after which the increment structure maps them back to \(\partial F/\partial X\) and \(\partial F/\partial Y\) [2509.10613]. The paper characterizes this dynamic-programming adjoint as exact because it differentiates the implemented operator itself.

## 4. Software architecture and integration with PyTorch

The library exposes a Python API backed by C++ and CUDA kernels, and all operations consume and produce PyTorch tensors [2509.10613]. Full compatibility with autograd means that both forward and backward passes participate directly in standard training loops [2509.10613]. This design aligns the library with typical deep learning workflows in which signature transforms or signature kernels appear inside a loss function or a trainable module.

The CPU code is described as exploiting SIMD and contiguous flattened memory layouts, while the GPU code uses shared-memory anti-diagonal rotation and block-of-32 tiling [2509.10613]. The library supports batched inputs and variable dyadic orders per input stream, \((\lambda_1,\lambda_2)\), to handle different sequence lengths efficiently [2509.10613]. It also performs lead–lag and time augmentation on-the-fly for signatures, which the paper states is faster and more memory-efficient than precomputing those transforms [2509.10613].

For high-dimensional paths, the computation of \(\Delta\) via `torch.bmm` is a key architectural choice, since it leverages batched matrix multiplication in PyTorch [2509.10613]. Likewise, the decision to avoid precomputing dyadically refined grids and instead reference refinement on-the-fly is intended to reduce memory use and improve throughput [2509.10613]. Numerical types follow PyTorch defaults; `float32` is identified as typical, and determinism is governed by standard floating-point behavior, with slight CPU/GPU differences possible due to parallel reduction order [2509.10613].

The paper includes minimal examples in which `psl.signature` computes flattened signatures for batched inputs, and `psl.signature_kernel` computes batched kernel values on CPU or GPU [2509.10613]. It also presents an MMD-style loss using the signature kernel:
\[
\mathrm{MMD}^2 = \mathbb{E}[k(R,R)] + \mathbb{E}[k(G,G)] - 2\mathbb{E}[k(R,G)],
\]
followed by `backward()` so that gradients flow to generated paths through the exact differentiation scheme [2509.10613]. A plausible implication is that the library is intended not only for feature extraction but also for direct use as an optimisation primitive in generative modeling.

## 5. Complexity, scalability, and benchmarked performance

The stated asymptotic complexity for truncated signatures of degree \(M\) on a piecewise linear path with \(L\) segments in dimension \(d\) is \(O(Ld^M)\) for the direct method [2509.10613]. Horner’s rearrangement reduces the number of multiplications and memory accesses but preserves the same asymptotic dependence on \(d^M\) [2509.10613]. This confirms a standard limitation of explicit truncated signatures: the method becomes expensive as either degree or ambient dimension increases.

For signature kernels on sequences of lengths \(L_X\) and \(L_Y\) in dimension \(d\), runtime is said to be dominated by two terms: \(\Delta\) computation with cost \(O(dL_XL_Y)\), and the PDE sweep with cost \(O((2^{\lambda_1}L_X)(2^{\lambda_2}L_Y))\) [2509.10613]. Memory use includes storing \(\Delta\) as an \(L_X \times L_Y\) array; on GPU, the forward \(\hat K\) requires only three anti-diagonals, giving \(O(L_X+L_Y)\) working memory [2509.10613]. The paper notes that CPU grid storage can be similarly reduced or streamed.

The reported benchmarks were run on Windows 11 with an i7-13700H, 32GB RAM, RTX 4060, Python 3.9, CUDA 11.8, and pySigLib v0.2.0, using the minimum over 50 runs [2509.10613]. The results are summarised below.

| Task | Baseline | pySigLib |
|---|---:|---:|
| Truncated signatures, forward CPU serial, \((128, L=256, d=4, N=6)\) | esig 1.1310s; iisignature 0.4104s | 0.0482s |
| Truncated signatures, forward CPU serial, \((128, L=512, d=8, N=5)\) | esig 11.4814s; iisignature 4.7908s | 0.3673s |
| Truncated signatures, forward CPU serial, \((128, L=1024, d=16, N=4)\) | esig 34.7836s; iisignature 14.3296s | 1.1512s |
| Truncated signatures, forward CPU parallel, \((128,256,4,6)\) | signatory 0.0558s | 0.0110s |
| Truncated signatures, forward CPU parallel, \((128,512,8,5)\) | signatory 0.4512s | 0.0896s |
| Truncated signatures, forward CPU parallel, \((128,1024,16,4)\) | signatory 2.2121s | 0.2988s |
| Truncated signatures, backward CPU parallel, \((128,256,4,6)\) | signatory 0.5918s | 0.1212s |
| Truncated signatures, backward CPU parallel, \((128,512,8,5)\) | signatory 5.9214s | 1.6595s |
| Truncated signatures, backward CPU parallel, \((128,1024,16,4)\) | signatory 18.2688s | 5.5374s |
| Signature kernels forward CPU, dyadic order 0, \((128,256,8)\) | sigkernel 0.3610s | 0.0118s |
| Signature kernels forward CPU, dyadic order 0, \((128,512,16)\) | sigkernel 1.3608s | 0.0364s |
| Signature kernels forward CPU, dyadic order 0, \((128,1024,32)\) | sigkernel 6.0874s | 0.2325s |
| Signature kernels forward GPU, dyadic order 0, \((128,256,8)\) | sigkernel 0.0088s | 0.0034s |
| Signature kernels forward GPU, dyadic order 0, \((128,512,16)\) | sigkernel 0.0376s | 0.0117s |
| Signature kernels forward GPU, dyadic order 0, \((128,1024,32)\) | sigkernel failed due to thread/memory limits | 0.0653s |
| Signature kernels backward CPU, \((128,256,8)\) | 2.3248s | 0.0322s |
| Signature kernels backward CPU, \((128,512,16)\) | 15.0445s | 0.1276s |
| Signature kernels backward CPU, \((128,1024,32)\) | sigkernel failed | 0.6019s |
| Signature kernels backward GPU, \((128,256,8)\) | 0.1475s | 0.0057s |
| Signature kernels backward GPU, \((128,512,16)\) | 15.3843s | 0.0231s |
| Signature kernels backward GPU, \((128,1024,32)\) | sigkernel failed | 0.1112s |

The paper interprets these results as showing substantial speedups, often greater than \(10\times\) on CPU and \(2\times\) to \(20\times\) on GPU, together with robust scalability to longer sequences and larger batches [2509.10613]. It also states that exact gradients improve training reliability [2509.10613]. This suggests that pySigLib’s performance gains derive not from a single kernel optimization but from the combination of memory layout, on-the-fly refinement, efficient batched \(\Delta\) computation, and the adjoint scheme that avoids a second PDE solve.

## 6. Applications, limitations, and practical considerations

The paper places pySigLib within several application domains for signature kernels, including hypothesis testing, distribution regression, path-dependent PDE solvers, and infinite-width kernel limits for deep sequence models [2509.10613]. In generative modeling for time series, especially in quantitative finance, signature kernels are described as powerful losses, including MMD and kernel scores, because they capture rich path structure [2509.10613]. pySigLib’s acceleration and exact gradients are said to make large batches and long sequences tractable, while improving reliability because backward passes no longer depend on approximate PDE gradients [2509.10613].

The paper also emphasizes practical preprocessing choices. In financial contexts, on-the-fly lead–lag and time augmentation are reported to provide substantial speedups and simplify preprocessing [2509.10613]. If sampling is highly irregular, time augmentation that includes timestamps can restore timing information in the features [2509.10613]. These points are implementation-oriented rather than theoretical, but they matter for empirical deployment.

The stated limitations follow directly from the underlying computational structure. Truncated signatures scale roughly with \(d^M\), so high degrees and large ambient dimensions can become prohibitive; the paper recommends signature kernels for high-dimensional tasks [2509.10613]. Higher dyadic orders \((\lambda_1,\lambda_2)\) improve PDE resolution but increase compute and memory, and the suggested practical strategy is to start with \(\lambda=0\) and increase only as needed [2509.10613]. Very large truncation degrees or aggressive kernel weightings such as large \(\lambda\) in \(\Lambda_m=\lambda^m/m!\) can cause growth in norms, and normalizing inputs, standardizing increments, or using time augmentation are listed as possible mitigations [2509.10613].

Parameter selection is framed as a balance between expressivity and cost. For signatures, the paper states that \(M=3\)–\(6\) is often practical; for kernels, it advises relying on PDE discretization and dyadic order rather than explicit \(M\) [2509.10613]. Device choice also depends on workload size: GPU use is recommended when batches are sufficiently large to saturate the block-of-32 tiling, whereas CPU may remain competitive on small workloads [2509.10613]. These remarks are not universal prescriptions, but they delineate the operational envelope within which the library was designed to perform well.

## 7. Relation to existing software and availability

pySigLib is explicitly compared with several existing libraries in the paper [2509.10613]. Against esig and iisignature, it is described as substantially faster for forward and backward signature computations, with the additional observation that iisignature recomputes signatures in the backward pass, leading to substantially longer runtimes [2509.10613]. Against Signatory, which provides GPU-enabled signatures and log-signatures with Horner-style updates, pySigLib’s CPU parallel and memory optimizations are reported to yield consistent speedups, and backward passes are faster in the reported benchmarks [2509.10613].

For signature-kernel toolkits such as sigkernel and sigkerax, the paper identifies two limitations of earlier approaches: precomputation of refined grids and constraints imposed by GPU thread counts or CPU memory [2509.10613]. pySigLib’s anti-diagonal shared-memory scheme and block-of-32 tiling are presented as removing thread-count limitations and scaling to longer sequences, while the exact differentiation scheme is said to outperform approximate PDE-gradient methods in both speed and accuracy [2509.10613]. The paper also distinguishes pySigLib from torchcde and torchsig, which it describes as focused on controlled differential equations and signature-inspired models rather than direct high-performance computation of signature transforms and signature kernels [2509.10613].

The paper summarizes the library’s distinctive features as GPU anti-diagonal tiling with shared memory, batched matrix multiplication for \(\Delta\), on-the-fly dyadic refinement and path transforms, and an exact, efficient gradient scheme obtained by differentiating the solver itself [2509.10613]. pySigLib is reported as open source at `https://github.com/daniil-shmelev/pySigLib`, installable via `pip`, with documentation at `https://pysiglib.readthedocs.io` [2509.10613]. Reproducibility information includes the reported software and hardware versions, and the paper notes that deterministic tensor operations are available subject to standard floating-point and parallel-reduction effects [2509.10613].

Source: https://www.emergentmind.com/topics/pysiglib