pySigLib: Fast Signature Computation

• pySigLib is a Python library for signature-based computation on sequential data, offering efficient truncated signatures and signature kernels.
• It leverages advanced CPU and GPU algorithms including on-the-fly dyadic refinement, SIMD-friendly memory layouts, and anti-diagonal tiling for improved throughput.
• The library integrates with PyTorch’s autograd to provide exact gradient computations, making it scalable for large-scale time series and generative modeling.

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” (Shmelev et al., 12 Sep 2025). 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 (Shmelev et al., 12 Sep 2025). 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 (Shmelev et al., 12 Sep 2025).

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 (Shmelev et al., 12 Sep 2025). 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 (Shmelev et al., 12 Sep 2025). 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 (Shmelev et al., 12 Sep 2025).

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 (Shmelev et al., 12 Sep 2025). 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 (Shmelev et al., 12 Sep 2025).

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 (Shmelev et al., 12 Sep 2025). 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$ (Shmelev et al., 12 Sep 2025). 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 (Shmelev et al., 12 Sep 2025).

The signature kernel between two paths $X : [0,T] \to \mathbb{R}^d$0 and $X : [0,T] \to \mathbb{R}^d$1 is the inner product in the tensor algebra. With weights $X : [0,T] \to \mathbb{R}^d$2, such as exponential weighting $X : [0,T] \to \mathbb{R}^d$3, the truncated kernel is

$X : [0,T] \to \mathbb{R}^d$4

where the inner product is the natural one on tensor powers (Shmelev et al., 12 Sep 2025). The paper emphasizes that pySigLib does not compute $X : [0,T] \to \mathbb{R}^d$5 and $X : [0,T] \to \mathbb{R}^d$6 explicitly to obtain this kernel. Instead, it uses the Goursat PDE characterization

$X : [0,T] \to \mathbb{R}^d$7

with boundary conditions $X : [0,T] \to \mathbb{R}^d$8, and discretizes this PDE on dyadically refined grids (Shmelev et al., 12 Sep 2025).

For a practical discretization over grid indices $X : [0,T] \to \mathbb{R}^d$9, with increments

$M$0

the finite-difference update is

$M$1

where

$M$2

Dyadic refinement orders $M$3 define a grid $M$4 of size $M$5 (Shmelev et al., 12 Sep 2025). 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 $M$6 in one flattened, contiguous array and updates each level in reverse order:

$M$7

Reverse-order updates permit in-place writes without temporary buffers, which minimizes allocations and memory traffic (Shmelev et al., 12 Sep 2025).

The second implementation is Horner’s method, described as signatory-style. The update is rearranged as

$M$8

The implementation separates an intermediate

$M$9

and then performs $S(X) = (S^0(X), S^1(X), \ldots, S^M(X)),$0 in-place, reusing a single preallocated contiguous buffer sized for $S(X) = (S^0(X), S^1(X), \ldots, S^M(X)),$1 and writing products in reverse order so that overwritten values are no longer needed (Shmelev et al., 12 Sep 2025). 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 (Shmelev et al., 12 Sep 2025).

For signature kernels, the CPU implementation first precomputes

$S(X) = (S^0(X), S^1(X), \ldots, S^M(X)),$2

and similarly $S(X) = (S^0(X), S^1(X), \ldots, S^M(X)),$3, then computes

$S(X) = (S^0(X), S^1(X), \ldots, S^M(X)),$4

via batched matrix multiplication using torch.bmm; the paper states that this step dominates runtime for large $S(X) = (S^0(X), S^1(X), \ldots, S^M(X)),$5 (Shmelev et al., 12 Sep 2025). Dyadic refinement is then applied on-the-fly by referencing $S(X) = (S^0(X), S^1(X), \ldots, S^M(X)),$6 at each stencil step, rather than precomputing refined grids (Shmelev et al., 12 Sep 2025).

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 (Shmelev et al., 12 Sep 2025). 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 (Shmelev et al., 12 Sep 2025).

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 (Shmelev et al., 12 Sep 2025). pySigLib instead differentiates the finite-difference solver itself to produce exact gradients in one backward pass over the same grid, reusing forward states (Shmelev et al., 12 Sep 2025). For $S(X) = (S^0(X), S^1(X), \ldots, S^M(X)),$7, the adjoint recurrences are given by

$S(X) = (S^0(X), S^1(X), \ldots, S^M(X)),$8

and

$S(X) = (S^0(X), S^1(X), \ldots, S^M(X)),$9

For dyadic refinement, multiple $$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},$$0 entries may depend on a given $$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},$$1, and the corresponding contributions are accumulated (Shmelev et al., 12 Sep 2025). Gradients with respect to the increments satisfy

$$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},$$2

after which the increment structure maps them back to $$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},$$3 and $$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},$$4 (Shmelev et al., 12 Sep 2025). 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 (Shmelev et al., 12 Sep 2025). Full compatibility with autograd means that both forward and backward passes participate directly in standard training loops (Shmelev et al., 12 Sep 2025). 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 (Shmelev et al., 12 Sep 2025). The library supports batched inputs and variable dyadic orders per input stream, $$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},$$5, to handle different sequence lengths efficiently (Shmelev et al., 12 Sep 2025). 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 (Shmelev et al., 12 Sep 2025).

For high-dimensional paths, the computation of $$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},$$6 via torch.bmm is a key architectural choice, since it leverages batched matrix multiplication in PyTorch (Shmelev et al., 12 Sep 2025). 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 (Shmelev et al., 12 Sep 2025). 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 (Shmelev et al., 12 Sep 2025).

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 (Shmelev et al., 12 Sep 2025). It also presents an MMD-style loss using the signature kernel: $$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},$$7 followed by backward() so that gradients flow to generated paths through the exact differentiation scheme (Shmelev et al., 12 Sep 2025). 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 $$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},$$8 on a piecewise linear path 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},$$9 segments in dimension $S^0(X)=1$0 is $S^0(X)=1$1 for the direct method (Shmelev et al., 12 Sep 2025). Horner’s rearrangement reduces the number of multiplications and memory accesses but preserves the same asymptotic dependence on $S^0(X)=1$2 (Shmelev et al., 12 Sep 2025). 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 $S^0(X)=1$3 and $S^0(X)=1$4 in dimension $S^0(X)=1$5, runtime is said to be dominated by two terms: $S^0(X)=1$6 computation with cost $S^0(X)=1$7, and the PDE sweep with cost $S^0(X)=1$8 (Shmelev et al., 12 Sep 2025). Memory use includes storing $S^0(X)=1$9 as an $X_{t_i}$0 array; on GPU, the forward $X_{t_i}$1 requires only three anti-diagonals, giving $X_{t_i}$2 working memory (Shmelev et al., 12 Sep 2025). 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 (Shmelev et al., 12 Sep 2025). The results are summarised below.

Task	Baseline	pySigLib
Truncated signatures, forward CPU serial, $X_{t_i}$3	esig 1.1310s; iisignature 0.4104s	0.0482s
Truncated signatures, forward CPU serial, $X_{t_i}$4	esig 11.4814s; iisignature 4.7908s	0.3673s
Truncated signatures, forward CPU serial, $X_{t_i}$5	esig 34.7836s; iisignature 14.3296s	1.1512s
Truncated signatures, forward CPU parallel, $X_{t_i}$6	signatory 0.0558s	0.0110s
Truncated signatures, forward CPU parallel, $X_{t_i}$7	signatory 0.4512s	0.0896s
Truncated signatures, forward CPU parallel, $X_{t_i}$8	signatory 2.2121s	0.2988s
Truncated signatures, backward CPU parallel, $X_{t_i}$9	signatory 0.5918s	0.1212s
Truncated signatures, backward CPU parallel, $z = X_b - X_a$0	signatory 5.9214s	1.6595s
Truncated signatures, backward CPU parallel, $z = X_b - X_a$1	signatory 18.2688s	5.5374s
Signature kernels forward CPU, dyadic order 0, $z = X_b - X_a$2	sigkernel 0.3610s	0.0118s
Signature kernels forward CPU, dyadic order 0, $z = X_b - X_a$3	sigkernel 1.3608s	0.0364s
Signature kernels forward CPU, dyadic order 0, $z = X_b - X_a$4	sigkernel 6.0874s	0.2325s
Signature kernels forward GPU, dyadic order 0, $z = X_b - X_a$5	sigkernel 0.0088s	0.0034s
Signature kernels forward GPU, dyadic order 0, $z = X_b - X_a$6	sigkernel 0.0376s	0.0117s
Signature kernels forward GPU, dyadic order 0, $z = X_b - X_a$7	sigkernel failed due to thread/memory limits	0.0653s
Signature kernels backward CPU, $z = X_b - X_a$8	2.3248s	0.0322s
Signature kernels backward CPU, $z = X_b - X_a$9	15.0445s	0.1276s
Signature kernels backward CPU, $$S(\text{segment}) = \sum_{k=0}^{\infty} \frac{z^{\otimes k}}{k!},$$0	sigkernel failed	0.6019s
Signature kernels backward GPU, $$S(\text{segment}) = \sum_{k=0}^{\infty} \frac{z^{\otimes k}}{k!},$$1	0.1475s	0.0057s
Signature kernels backward GPU, $$S(\text{segment}) = \sum_{k=0}^{\infty} \frac{z^{\otimes k}}{k!},$$2	15.3843s	0.0231s
Signature kernels backward GPU, $$S(\text{segment}) = \sum_{k=0}^{\infty} \frac{z^{\otimes k}}{k!},$$3	sigkernel failed	0.1112s

The paper interprets these results as showing substantial speedups, often greater than $$S(\text{segment}) = \sum_{k=0}^{\infty} \frac{z^{\otimes k}}{k!},$$4 on CPU and $$S(\text{segment}) = \sum_{k=0}^{\infty} \frac{z^{\otimes k}}{k!},$$5 to $$S(\text{segment}) = \sum_{k=0}^{\infty} \frac{z^{\otimes k}}{k!},$$6 on GPU, together with robust scalability to longer sequences and larger batches (Shmelev et al., 12 Sep 2025). It also states that exact gradients improve training reliability (Shmelev et al., 12 Sep 2025). 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 $$S(\text{segment}) = \sum_{k=0}^{\infty} \frac{z^{\otimes k}}{k!},$$7 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 (Shmelev et al., 12 Sep 2025). 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 (Shmelev et al., 12 Sep 2025). 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 (Shmelev et al., 12 Sep 2025).

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 (Shmelev et al., 12 Sep 2025). If sampling is highly irregular, time augmentation that includes timestamps can restore timing information in the features (Shmelev et al., 12 Sep 2025). 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 $$S(\text{segment}) = \sum_{k=0}^{\infty} \frac{z^{\otimes k}}{k!},$$8, so high degrees and large ambient dimensions can become prohibitive; the paper recommends signature kernels for high-dimensional tasks (Shmelev et al., 12 Sep 2025). Higher dyadic orders $$S(\text{segment}) = \sum_{k=0}^{\infty} \frac{z^{\otimes k}}{k!},$$9 improve PDE resolution but increase compute and memory, and the suggested practical strategy is to start with $$S(x * y)_{[a,c]} = S(x)_{[a,b]} \otimes S(y)_{[b,c]}.$$0 and increase only as needed (Shmelev et al., 12 Sep 2025). Very large truncation degrees or aggressive kernel weightings such as large $$S(x * y)_{[a,c]} = S(x)_{[a,b]} \otimes S(y)_{[b,c]}.$$1 in $$S(x * y)_{[a,c]} = S(x)_{[a,b]} \otimes S(y)_{[b,c]}.$$2 can cause growth in norms, and normalizing inputs, standardizing increments, or using time augmentation are listed as possible mitigations (Shmelev et al., 12 Sep 2025).

Parameter selection is framed as a balance between expressivity and cost. For signatures, the paper states that $$S(x * y)_{[a,c]} = S(x)_{[a,b]} \otimes S(y)_{[b,c]}.$$3–$$S(x * y)_{[a,c]} = S(x)_{[a,b]} \otimes S(y)_{[b,c]}.$$4 is often practical; for kernels, it advises relying on PDE discretization and dyadic order rather than explicit $$S(x * y)_{[a,c]} = S(x)_{[a,b]} \otimes S(y)_{[b,c]}.$$5 (Shmelev et al., 12 Sep 2025). 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 (Shmelev et al., 12 Sep 2025). 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 (Shmelev et al., 12 Sep 2025). 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 (Shmelev et al., 12 Sep 2025). 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 (Shmelev et al., 12 Sep 2025).

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 (Shmelev et al., 12 Sep 2025). 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 (Shmelev et al., 12 Sep 2025). 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 (Shmelev et al., 12 Sep 2025).

The paper summarizes the library’s distinctive features as GPU anti-diagonal tiling with shared memory, batched matrix multiplication for $$S(x * y)_{[a,c]} = S(x)_{[a,b]} \otimes S(y)_{[b,c]}.$$6, on-the-fly dyadic refinement and path transforms, and an exact, efficient gradient scheme obtained by differentiating the solver itself (Shmelev et al., 12 Sep 2025). pySigLib is reported as open source at https://github.com/daniil-shmelev/pySigLib, installable via pip, with documentation at https://pysiglib.readthedocs.io (Shmelev et al., 12 Sep 2025). 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 (Shmelev et al., 12 Sep 2025).