Scalable Signature Kernel Computations via Local Neumann Series Expansions
Authors: Matthew Tamayo-Rios
Alexander Schell
Rima Alaifari
Venue: NeurIPS 2025 poster ()
Keywords: signature kernel
goursat
rough path
machine learning
time series
Primary area: general_machine_learning
TL;DR: A new scalable accurate approach to computing signature kernels using Neumann-series expansions

Abstract: The signature kernel is a recent state-of-the-art tool for analyzing high-dimensional sequential data, valued for its theoretical guarantees and strong empirical performance. In this paper, we present a novel method for efficiently computing the signature kernel of long, high-dimensional time series via adaptively truncated recursive local power series expansions. Building on the characterization of the signature kernel as the solution of a Goursat PDE, our approach employs tilewise Neumann‐series expansions to derive rapidly converging power series approximations of the signature kernel that are locally defined on subdomains and propagated iteratively across the entire domain of the Goursat solution by exploiting the geometry of the time series. Algorithmically, this involves solving a system of interdependent Goursat PDEs via adaptively truncated local power series expansions and recursive propagation of boundary conditions along a directed graph in a topological ordering. This method strikes an effective balance between computational cost and accuracy, achieving substantial performance improvements over state-of-the-art approaches for computing the signature kernel. It offers (a) adjustable and superior accuracy, even for time series with very high roughness; (b) drastically reduced memory requirements; and (c) scalability to efficiently handle very long time series (e.g., with up to one million data points or more) on a single GPU. As demonstrated in our benchmarks, these advantages make our method particularly well-suited for rough-path-assisted machine learning, financial modeling, and signal processing applications involving very long and highly volatile sequential data.

Scalable Signature Kernel Computations via Local Neumann Series Expansions
Matthew Tamayo-Rios
Alexander Schell
Seminar for Applied Mathematics & Al Center ETH Zurich, Switzerland
Department of Mathematics
Technical University of Munich, Germany
Rima Alaifari
Department of Mathematics
RWTH Aachen University, Germany
Abstract
The signature kernel [10 is recent state-of-the-art tool for analyzing high- dimensional sequential data, valued for its theoretical guarantees and strong empiri- cal performance. this paper, we present novel method for efficiently computing the signature kernel of long, high sigh-dimenssonal time series via adaptively truncated recursive local power series expansions. Building on the characterization of the signature kernel as the solution of Goursat PDE [17]. our approach employs tilewise Neumann-series expansions approximations of the signature kernel that are locally defined on subdomains and propagated iteratively across the entire domain of the Goursat solution by exploiting the geometry of the time series. Algorithmically, this involves solving system of interdependent Goursat PDEs via adaptively truncated local power series expansions and recursive propagation of boundary conditions along directed graph in a topological ordering. This method strikes an effective balance between computational cost and accuracy, achieving substantial performance improvements over state-of-the-ant art approaches for computing the signature kernel. It offers (a) ad- justable and superior accuracy, even for time series with very high roughness; (b) drastically reduced memory requirements; and (c) scalability to efficiently han- dle very long time series (one million data points more) demonstrated our benchmarks, these advantages make our method particularly well-suited for nough-path-assitted machine learning, financial modeling, and sig- nal processing applications involving very long and highly volatile sequential data.
derive rapidly converging power series
single GPU. As
Introduction
Time series data is ubiquitous contemporary data science and machine learning, appearing in diverse applications such satellite communication, radio astronomy, health monitoring, climate analysis, and language or video processing, among many others (21]. The sequential nature of this data presents unique challenges, it is characterised by temporal dependencies and resulting structural patterns that must be captured efficiently to model and predict time-dependent systems and phenomena with accuracy. Robust and scalable tools for handling such data in their full temporal complexity are thus essential for advancing machine learning applications across these domains. One particularly powerful approach in this direction is the signature kernel [TO] [7]. the Gram matrix of high-fidelity feature embedding rooted in rough path theory and stochastic analysis 14 g). which has gained relevance as an increasingly popular tool in the modern analysis of sequential data
13
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

Conceptually, the signature kernel (of family of multidimensional time series) is the Gram matrix of the signature transform, a highly informative and faithful feature map that embeds time series into a Hilbert space via their iterated integrals, thus uniquely capturing their geometry and essential time- global characteristics in hierarchical manner (9]. This structure situates the signature kernel naturally within the framework of reproducing kernel Hilbert spaces (RKHS) and enables its practical utility through kernel learning techniques. The rich intrinsic properties of the underlying signature transform further confer several strong theoretical properties to the signature kernel, including invariance under irregular sampling and reparametrization, universality and characteristicness, and robustness to noise [6.. Its wide-ranging theoretical interpretability combined with its strong weal-world efficiency have elevated the signature kernel to a state-of-the-al tool for the analysis of dependeenndnt data ([3].
However, existing methods for computing the signature kernel suffer from significant scalability issues, particularly when dealing with long or highly variable time series ! The reason is that these methods typically either approximate truncations of the signature transforms via dynamic programming [10] or solve global Goursat PDE for the signature kernel using two-dimensional FDM discretizations [77, which both involve at least quadratic storage complexity relative to the length E of the time series, thus resulting in prohibitive memory usage for long or very rough time series. These limitations, and the quadratic complexity in the number of constituent time series (a common kernel method bottleneck not specific to the signature kernel,, can severely obstruct the application of signature kernels to large-scale, real-workee datasets; see, e.g.. [19 and the references therein.
To address this, we introduce novel approach to compute the signature kernel based on adaptively truncated local (tilewise) Neumann series expansions of the solution to its characterizing Goursat PDE; see Figure 1 for an overview. Partitioning the domain [0, 1]2 of this PDE solution into tiles according to the data points of the input time series, we derive rapidly convergent power series representations of the solution each tile, with boundary data obtained from adjacent tiles. This enables efficient computation of the signature kernel in terms of memory and runtime, while ensuring superior accuracy and significantly improved scalability with respect to the length and roughness of the kernel's 's constituent time series. Since only local expansions are wered-rather than a global & x E grid-memory usage grows much more slowly with E, allowing the method to handle very long time series, even up to a million points each on a single GPU, where PDE or solvers typically run out of memory. By explicitly leveraging the time-domain geometry of piecewise- linear time series in this way, our method supports localized computations of the signature kernel that both highly parallelizable and memory memory-efficient, even for very long and rough time series.
progriimnggg----g----.. based
Specifically, our main contributions are:
We propose a tilewise integral-equation approach to construct recursive local power series expansions of the signature kernel, based on Souudday-ttoboundde leverages the kernel's time-domain geometry to enable its computation with significantly reduced memory requirements and computational cost without compromising accuracy.
Neumann Series for the Signature Kernel:
propagation and paired with an adaptive truncation strategy. This
2. Parallelizable Local Computation: Exploiting the piecewise-linear structure of the input time series, we partition the full kernel domain into ordered tiles supporting parallel local Neumann expansions with adjustable precision and minimal global communication.
3. Scalability to Very Long and Rough Time Series: Our method, termed PowerSig, achieves scalability to very long (over 10th points on a single GPU) and highly volatile (rough) time series, addressing key limitations of existing methods.
4. Empirical Validation On Benchmarks: We demonstrate the practical advantages of our method through comprehensive benchmarks against several state-of-thr ofttt--arttt signature kernel solvers, demonstrating superior accuracy, runtime, and memory efficiency of our approach.
The remainder of this paper elaborates on these contributions, starting with a brief review of re- lated work (Section I.I). followed by a detailed description of our methodology (Section 2 and a presentation of our experimental results (Section 3. Proofs are given in Appendix A.
'Variability of a (discrete-time) time series is quantified by the sum of squared differences between consecu- tive points of the time series and referred to as its 'roughness'.

---

The codebase for our method, including all implementation details, is provided in the supplementary material and publicly available at: hnttps://gita .tt/:///aebbbbbs////rrrsss....i.
1.1 Related Work
The computation of Gram matrices for genature-transformee time series was first systematically studied by Király and Oberhauser [10), who identified the signature kernel as foundational link between rough path theory, data science, and machine learning an interplay broadly envisioned earlier by Lyons [[2]. Chevyrev and Oberhauser [6] further extended these theoretical foundations, and introduced a statistically robust variant of the signature kernel via appropriate time series scaling.
Computationally, Salvi al. [17] significantly advanced beyond earlier dynamic programming approaches for computing truncated kernels [0]] by characterizing the (untruncated) signature kernel through linear second-order hyperbolic (Goursat) PDE. This approach, implemented in widely-used libraries such as sig sigkernel and GPU-accelerated KSig [19], provided efficient signature kernel computation through diniteedifference PDE approximations.
However, PDE-based methods, despite providing highly parallelizable and fast routines for short- to moderate-length time series, exhibit poor scalability due quadratic memory usage, becoming impractical for long or rough time series beyond a few thousand to tens of thousands of time steps
[ (our experiments confirm computational limits of approximately 10s steps for sigkernel and 16 x 10 for KSig on 4090 RTX GPUs). Dynamic programming remains an option for computing kernels of truncated signatures 110, although it likewise suffers from poor scalability. Recent efforts to reduce this cost through random Fourier features or other low-rank approximations (20) often degrade in accuracy for larger time-series length or time series of higher dimension.
19
Our proposed method, PowerSig, circumvents these issues through tilewise local power series expansions of the signature kernel, avoiding the need for global storage footprint that scales quadratically with time series length. By viewing the Goursat PDE as a Volterra integral equation, compute rapidly convergent Neumann series expansions of the kemel locally, storing only series coefficients rather than full two-dimensional arrays. This localized strategy enables fast and accurate signature kernel computations for extremely long (over 10 points) and high-dimensional time series on single GPUs, significantly improving scalability and efficiency.
During the preparation of this manuscript, came across a preprint by Cass al. (5] which appeared concurrently to the release of our first version. Their work is also based on recursive power series expansions of the signature kernel, albeit different mathematical and algorithmic premises. An empirical and conceptual comparison between our method and theirs is provided in Section B.2
2 Signature Kernels via Recursive Local Neumann Series
This section presents our tilewise Neumann-series expansion method for the signature kernel. We begin by revisiting the PDE characterization of the signature kernel [17 (Section 2.1 and recast it as an equivalent integral equation that allows for a recursive, tile-based decomposition of its solution with only boundary values exchanged between adjacent tiles (Section 2.2). This local formulation produces rapidly convergent nested power series expansions on each tile and supports an adaptive truncation scheme to balance aceuracy and computational cost (Section 2.3). By storing and passing only local series coefficients, the method achieves substantial memory savings compared to global PDE or dynamic- programming solvers (Section 3). Figure summarizes the core idea.
2.1 The Signature Kernel
Let æ = (æ.,..., 23) ( C IRd and y (y),..,,, yr) C Rd be time series of common length E E N (if lengths differ, one may pad the shorter time series by repeating its final entry). For any such time series = (Z1, zz), define its affine interpolami 2 : [0, 1] + Rd by
2(t) = + (e- 1) (t (t- = AAz, tE - [
11,, 1),
where A&* = x (This is just the unique continuous piecewise linear function interpolating the points 21,..., 24.) The derivative of Σ (defined almost everywhere on [0, 1)) given by
2'(1) = (f - 1) A2z, LE 8-1 (
3

---

For later use, we now partition the unit square [0, 1|2 into riles as follows. Set
and 1
and define the open tiles and their closures (in the Euclidean topology on IR22) as
P=k==(Ok- xx+1) x (7i, 7+1) and Tk = Da,1 = (o x++1 1| x [77 7+1l
(1)
The tiled domain is then given by D = - DR., with boundary 0D - [0, 1] D.
Finally, let (.--) denote the Euclidean inner product on R Building the PDE characterization in our main object of interest as follows.
17 Theorem 2.5], we
Definition 2.1 (Signature Kernel). The signature kernel of and is the unique continuous function
K Kæ.y : [0, 1 + solving the hyperbolic (Goursat) boundary value problem
0 (s,
Oset. Pr,y (8,
(2)
K(O,-) ((,,0) 1,
where the coefficient function Pa.3 D + R is defined tilewise by Pa.y ( 1) = (ø' (s), / (1))
(s,t) Da.-
An equivalent formulation defines Pr.y on all of [0, 1|2 via
Pa.y [0, (s.t) + (e - 1)² Pk. (s, Pk. Ai Aky. 1-1=
(3)
In integral form, the boundary value problem is then equivalent to the Volterra integral equation
K(s,t) = So Pa (u, (u,v) K(u, v) du de, (s, E [0, 1J².
(4)
Standard fixed- point arguments (Appendix A.2 guarantee that 44 has a unique solution in C((O, 1)3). the space of continuous functions on [0, 1 This yields the following well-know result:
Proposition 2.2. The Goursat problem 2 has unique solution in C((0, 1)3); particular, the signature kernel of and y is well-defined.
The advantage of the Volterra formulation (4 is that it naturally leads recursive power series expansion for the signature kernel K, as we will now explain. A key step towards this observation, which is also used in the proof of Proposition 2.2 is the following lemma:
Lemma 2.3, For any bounded, measurable function @ + R defined on a closed rectangle == [ai, bi] (a2,b2 [0, 102 the integral operator Te (C(I), 11-1)al-1100 11x0) + (C(I), given by
(sf(s, = S
[ e(u, v) f(u, v) dudu, (s,t)
has spectral radius zero, that is, r(T,) supxeatre) [V] contimuous functions on 1, equipped with the supremum norm || == suP( [f(s, f)}-
Here, C(I) denotes the space of
This fact will justify the use of Neumann series when inverting operators of the form id Tp below.
2.2 A Recursive Local Power Series Expansion of the Signature Kernel
Our computational strategy is approximate the signature kernel K by a directed family of local power series expansions constructed recursively over the tiles (T.. In the spirit of the Adomian Decomposition Method (ADM) [ 22 we seek establish a representation
K(s,t) = Σ (s,))
(5)
'Introducing the signature kernel of and y as the solution to its characterizing PDE circumvents the need specify it as the inner product of the signature transforms of and y. originally formulated in [7].

---

Σ
Ey be-1zlo ) Aas-slin) n--u-1lon))
nod
T2.1
Na
Pec
75.5
75.
T2.1
13.1
Rocantre Look.
Adcettn Trika Thaeaaiie
Figure 1: Summary of our method ( (PowerSig) for computing the signature kernel of two time series via recursive local Neumann expansions. Panel A: The PDE 8 Pu.y induces (xr, d)-dennndnnt partition of [0, 1|2 into tiles Tk.1- Arrows indicate the sequential propagation of boundary conditions across tiles, with decreasing colour intensity corresponding to later propa- gation steps. Tiles receiving arrows of the same colour form groups whose local series be computed in parallel. Panel B: On each given tile Tk.1, the kernel admits the recursive local Neu- mann series expansion: ku, - TV, TW,,
K (ok 7:)
(65-r) Eij= (s (t-) Ti)², which converges uniformly on (s,t) E Tx,1- These tilewise expansions depend on boundary values from neighbouring tiles (Kk-1.1 and Kg,1-1), with arrows indicating the directions of integration from boundaries to the tile interior. Panel C illustrates adaptive series truncation and final kernel assembly. Tile shading intensity encodes local truncation depth, which is adaptively determined by the magnitude of Pk.1 (Øx+1 Ixx) (3+ y:). Darker tiles indicate the necessity for deeper (higher-ordss oorde) expansions, while lighter tiles allow shallower truncation.
with terms Kij 10, 1]2 R. that are easy compute and, for our purposes, take the form
(s, = Σ (s,
(s,
ck. (1.5) (s ox) 7))
(6)
k,1=1
with some tile-dependent coefficient sequences
Ck.1 = (c(t)) (i,5) E N3) E (NO),
(7)
which are summable and (for (k; * (1, 1)) defined recursively by Ck.1 = Φx. (c C--1.1) for maps Φx,1 l (N3) x2 E, (142) where 0 (resp. (1.1) depends only on its second (resp. first) argument.
(8)
The sets Te. are half-open, Tx.1 = [ Tk41) [T.. 7+1) for k. < é 1, and closed on the last row/column (i.e. when = - or 1 = - 1), TR.E-1 == Tke-1, P.-1, = T2-1. Thus (Tr.))e defines a partition of the Goursat domain [0, 1j2.
8-1
We organize the recursion K over all (e- 1) tiles, computing the coefficients 7 on each tile Te, via Neumann expansion from boundary data the adjacent tiles TR11;1 Tk,1-1 (whose coefficients (ck-1,1 and C2,1-1) are already available). Section 2.2.2 details the procedure, beginning with the
bottom-left tile Ti, in Section 2.2.1
The goals of this scheme are twofold: (a) choose the coefficients 7 so that the series 6 converge rapidly on each tile, giving power series representation of the signature kernel localizations
K = K Te. Jo 0b+1l X [T. 7441] 3 (s, t) - K ((s,t) E R.
and (b) truncate these tilewise series expansions 9 so as to obtaín numerically stable and efficient global approximation of the whole kernel K; see Section 2.3 for both.

---

2.2.1 Rapidly Convergent Power Series on the First Tile
On the first tile Ti,1 10, 02) [0, 72). we adopt 5 as an ansatz and assume
(1,1 - Σ (s,(s,t), レ,==0
where (i = cft,)) 5't
(10)
with Σis [(d)) < oo and only diagonal coefficients nonzero (c JTi Eis [i)((W) dw < co, and Fubini applied the integral equation gives
(i.)) for + j). Then
Kis with (s)(8,1) = Sol (u u) KK-15-- v) du Kig = Σ
(11)
+==0
i,j=0
and for Ko,0 = pointwise on 71, Since Pa.lTT.. = P1,1 (see 3 ), simple induction yields
pi
if = J. s'
(s, =
aay
pointwise on Ti.1-
if
+ J.
By the identity theorem for power series,
implies that the desired (41.1 (i.)) [ij No) must read
Pi.1 (2!)
11,1
(1.)
8
(12)
where dig is the Kronecker delta. On Ti, 1, the decomposition ansatz (T0 thus yields the well-known Lemma 2.4. On the first tile Ti,1, the signature kernel
has the form:
Pi. K(s,t) Σ s
{ Jo (2 (2Vleaals),
P1,1 <0,
uniformly in (8,t) E Ti, (13)
(it)
Lo (2.
P1,1 20.
where Jo and To are the Bessel and modified Bessel functions of the first kind of order 0, respectively.
The truncation error decays as (nn)-22) in the order n, making the series (3) highly effective for approximating K1,1, especially when |P1.1| I is moderate (larger |PL.1] need higher truncation orders).
2.2.2 Recursive Neumann Series for Propagating the Signature Kernel Across All Tiles
The recursion (8) for the local power series coefficients c₁2 (cl(.)) E (1())) the remaining tiles starts from the base coefficients (12) and proceeds as follows.
For k,
define the (propagation) operators Tk :CTk,) C(Tk.) by
S S'J"" Pu.if v) du (s, E Tsaa
(Tr.f)(s,t)- =
(14)
(cf. Lemma 2.3 and set To, == {0} x [m, T+ and Tk = lox, det1] x {0}-
Proposition 2.5. For each k 1, ...., 1, the restricted kernel KK.1 = Klr,, from 9 satisfies
(o..)) R-11 + (( - (2122 uniformly on Tx, KR.1 = Σ Th
(15)
n=0
(a.T) for Ki.l = (K,,(a,T) T) and the 'curried" functions Kig, Th.1 3 (u, + (o, v) and 15
Tx.1 3 (u,w) + ku,((u, Here, K0,1 = K K|TT, and K4 = K KlTTo and K0.0 = ((0,0) =
The identities in 15 yield the desired tilewise power-series representation coefficients 7 are determined recursively from those on the tiles immediately to the left and below. The following example illustrates this.
5 66 On each tile, the
This is assumption only for the moment will establish
10
provable identity in Lemma 2.4

---

Example 2.6 (Evaluating 15
on T1.1, T1.2. and T2,1). With
and No, = 14.0
No,0 = L, the recursion 15 gives on the bottom-left comerr tile Ti,1 that
K1,1
Tr
T7,1
T Σ T uniformly on Ti,1
(16)
n=0
where the operator Tí,1
((..1) PI, v) dudu) is applied repeatedly. Since
(s,,))(ss (o) for each E No (as can be readily verified by induction)) the recursion (16
(ntya
precisely recovers the expansion
13 On the adjacent tile Ti,, the identity 15) yields
80
(1,2(8,1)
Σ
T (0).0)
1,2
T2
(3.7a)
T 0,2
rx) (1,
PL.1 Th
AuPi (t J)li(e
T +
T Tr2 (s,
T
E
(it)²
(i
1)3
(.=0
where the last equality used the definition of and that, for each iteration index E No.
(s,,)*)(s = PL.2 []] TT (TT--1"" (*, u) du de =
Př,
inn T2)"
(17)
(i + 1)
with (z)* 11 i), as one verifies immediately (Lemma A.2). Analogous computations show
(2,1(s = Σ
(ke + l) - 1y 1y(s - x2
k.1=0
uniformly in (s,t)= 6 T2.1-
The observation 17 is recorded as Lemma A.2 for later use. Given the recursion (1). now need
an explicit algorithm (#) to extract the coefficients in (8) and thus build the approximations in (5).
The next section provides this (Propositions) 2.7 and 2.8
2.3 Computing the Neumann Series Coefficients
To algorithmically extract the power-series coefficients from tilewise Neumann recursions (5)).
we encode the action of the propagation operators Tu on the monomial basis {s*12 (i,j) E N C
((T..)) via a simple Vandermonde scheme: Define the power map 1) [0, t (No) by
((r) = (r* i No) (1, ...... and, for (cg)a520 E El (NG)), define the doubly-infing matrices C E (Lo (No), by
C= (c.)) a (az) (Ei20 :20
=.
(18)
By the Weierstrass M-test, each such C then induces continuous function Cla,)) [0, 1|2 3 (s, - (1(()) Cn(()) = Eij20 E,,
(19)
where ( denotes the dual pairing between too (No) and & (No). The localised kernels 9 can then all be represented in the explicit form (19) for some recursively related coefficient matrices C..l. Proposition 2.7. For each k, E {1,..., 1}. there is C.a. = (ck))..20 E (NG) such that
Kk,1 Caay (T.))
and ckl== -E in ((N)).
(20)
m=
with Kg.1 = limm- E Chi.l Isa.r
uniformly on Tw. The sequence (Ch.)) nzo C (N is
recursively defined by with initial value
Ch.. pr.iLoo, Chil Ch.lRr (n E No) Ohel (ayôio + B,000 - 780: dos)
(21)
(22)
7

---

for (o.)= Chil-1 (TT)) (Be) = (CK-Li)) n(ox) E (No) and q ( ((-)))((T)) and the boundary coefficients Cuu = C4. = (ooo--oo)))))) dos)).000000 each E No. Identities 21
use the matrices
Log = (T- and R = T(T G(r)),
1= (6.)).220 H(a) (026 G(T) = (T' ).520 s (
T (f:111
1,520'
Numerically, the dominant cost in applying (15) is the double integration S So i.c., the multiplica- tions Ro, in (TT) This cost depends on the expansion center (sk, 4i) in the representation
(8, E
(d st) - -4) ((x, E TR.))
(23)
Proposition 2.7 establishes (23 for (sk, (r) = (0, 0) and yields a direct implementation via (20 (2I). The next result shows that centering instead at the tile corner (sk, ft) (ox, 7)) sub- stantially reduces computation while preserving uniform convergence. In what follows, O de- notes the Hadamard product of doubly-infinite matrices (i.e., the entry-wise matrix multiplication (o1,5).,120 O (bi,) == (a., )2,520), (o.,bb,)).,20)). and we abbreviate E' =
Proposition 2.8. each k, E {1...., }.), the localised solution solutiio ( has the tile-centered expansion
[k,(8,1)= = Σ -
(1.0 4 (s (s-ox)) (x) TT) T for Ctk. ((cb)) (2(6.1)
1,500 = Ax.1 ④ BE, ② W, (24)
minnt.)) uniformly in (s,t) E Tris where W- = (Wi,))s20 with way = mix
well as
Ak. (Pk. min(i.))
1.200' and the matrix Bad = (b/) (b(b,)) 1220 is defined by (k:1)
0iitr = ar
(k.)) and 0,g-
(k.) == - (k.))
for each (i,r)
(25)
with sf1,1) = ar (1,1) = bu for each r E No. and recursively (recalling 18 for notation),
(a.h.)) r20 = (ih,(11)), m(1//e)))(1())))
and (B(h.)) r20 = [C nT(1/E).
(26)
Equations 24 (25). and (26) define an efficient, tile-centered implementation of (8. Proposition A.3 provides rigorous priori bounds for the induced Gram-matrix approximation error. The next section presents numerical evaluation of this method covering accuracy, memory usage, and runtime- and illustrates its applicability downstream tasks on real data.
3 Numerical Experiments
We evaluate our method, PowerSig, in terms of accuracy, memory usage, and runtime. Specifically, compute the self-signatur kernel of randomly drawn wwo-dimensional Brownian motion paths on [o, 1] at increasing sampling frequencies, using sample lengths & = 2 + for A 2 0, constrained only by GPU memory; comparisons are made against the state-of-the-art KSig library [19]. We use KSig both with its truncated signature kernel and with its PDE-based solver at the default dyadic order. Experiments were
NVIDIA RTX 4090 GPU (24 GB). Unless noted, PowerSig truncation order for the tile-center local series 24 is fixed at 7, although higher orders are equally feasible. Accuracy is reported as Mean Absolute Percentage Error (MAPE) relative the KSig truncated signature kernel (order 1, truncation level 21), memory is peak GPU usage, and runtime is total execution time. Each point averages 10 independent runs.
Accuracy. Figure 2 compares the accuracy between PowerSig and the PDE-based solver from KSig across wwo-dimensional Brownian motion paths of length up to 513 (the maximum length manageable by the truncated signature kernel; left panel) and on two-dimensional fractional Brownian motion paths of fixed length 51 across decreasing Hurst indices (from 0.4 down 0.005; right panel). PowerSig, despite employing only a modest truncation order, achieves superior accuracy and remarkably low error levels as length and irregularity ('roughness') of the input time series increase. PowerSig's g robust performance on time series with low Hurst indices further suggests significantly enhanced numerical stability of the method, particularly for highly irregular (rough) trajectories.
'Here and in the definition of I' we adopt the convention g 0.
8

---

10
10:
10'
10
1071
10'
35
Figure 2: Comparison of Mean Absolute Percentage Error (MAPE) between PowerSig and the PDE- based solver from KSig. Left: for two-dimensional Brownian motion paths [0, 1] across increasing path lengths Right: for two-dimensional fractional Brownian motions of fixed length 51 across increasingly irregular sample paths (decreasing Hurst index, swept through progressively rougher regimes); the right panel reports MAPE relative to the signature kernel truncated at level 180.
Memory Usage and Runtime. Figure 3 highlights the practical advantages PowerSig in terms of GPU memory usage and runtime. Specifically, its localized, tile-based computation drastically reduces memory overhead, enabling computations on paths of length = 524 289 with under 720 MB GPU memory, which is orders of magnitude lower than both PDE- and dynamic-progranme programming -baed methods. The inherent sparsity of the propagated Neumann-seriin expansions enables efficient mem- ory management, allowing PowerSig to handle substantially larger-scale problems and overcome the storage bottleneck associated with full-grid methods.
Mermayy Usage
12'
4230
3230
3230
Figure 3: Peak GPU memory usage (left) and runtime (right) for computing the signature kernel on two-dimensional Brownian motion paths, comparing PowerSig with the truncated-signatinn (KSig) and PDE-based (KSig PDE) solvers. PowerSig achieves substantially lower memory consumption and maintains computational feasibility for large well beyond the limits of the alternative methods.
Time Complexity PowerSig retains the O(Ed)) runtime scaling of existing methods but signif- icantly improves space complexity to O(EP), where P (polynomial truncation order) is typically much smaller than f. By storing only single diagonal of coefficient blocks, PowerSig enables the processing of much longer paths than previously feasible with existing approaches.
Empirical Evaluation Real and Large-Scale Settings
Beyond the above benchmarks, we also assess downstream performance and compare PowerSig with recent low-rank and random Fourier Feature (RFF) approximations. Unless noted otherwise, use the default truncation policy and report averages multiple independent runs. All figures,

---

implementation details, and full hyperparameter grids appear in the supplement (downstream figures have been moved to Section Bf the supplement to comply with page limits).
(A) Bitcoin price regression (Salvi el al. (I7I). Figure 4 shows train and test fits (two-day rolling average) for kernel-ridge regression (MAPE) on the public bitcoin pricing dataset featured in Salvi al. in [17]. On the test set, PowerSig attains 2.81% MAPE versus 3.23% for the (RBF-assisted) KSig-PDE. For the underlying Gram-matrix construction, peak memory for KSig-PDE scales ()) (with N windows and window length f), whereas PowerSig uses only ( For the present setup (N = 299, L = 36), this extrapolates to roughly 1.6,GB for KSig-PDE versus ~.0.038,MM for PowerSig. This
illustrates that the near-exts exact regime enabled by PowerSig's ' linear-in in-length memory profile yields tangible predictive gains at far lower resource cost, with particularly clear benefits even at short window lengths.
(B) UEA Eigenworms classification and RFF/lowrrak baselines. We benchmark PowerSig and KSig-PDS: against linear/RBF kernel SVMs and the recent specialized RFF-based method RFSF- TRP from [20] on the standard Eigenworms dataset with input window lengths L {16, 32....
... ... 1024]. As shown in Figure 5 PowerSig (and KSig-PDS up to L = 128 before OOM) remains competitive and rises to 61.1% accuracy at L = 1024, whereas RPSF-TRP attains a slightly higher peak of 62.5% at L 1 28 but exhausts memory for larger I., consistent with the storage advantages in Figure 3 These results show that narrow performance gaps often ascribed to inductive bias while maintaining feasibility.
substantially extending the input window enabled here at scale by PowerSig-
(C) Long-horizon periodic signals (industrial/sensing proxy). Motivated by predictive main- tenance (near-periodic gearbox/turbine vibrations) and narrow-band I/Q radio signals, we generate synthetic near-periodic time series with adjustable period length. For representative instances (see, e.g., Figure (7). SVM- reggesssssionn error decreases monotonically as input windows span multiple periods. As shown in Figure(6] PowerSig sustains this behavior at window lengths beyond the reach of conventional or low-rank gegnature-kernd methods, while peak memory grows linearly with window length.
(D) High-dimensional scaling. Complementing Figure we fix path length = 4096 and vary dimension d from 2 to 8192. (JsoowsSswwws stable accuracy and pear-perfect linear runtime from d = 64 to 4096 (slightly sublinear outside), with memory following our one-strip tiling profile. This corroborates the practicality of PowerSig in high-dimensional sensing and multivariate finance.
Overall, across real regression and classification tasks and stress tests in length and dimension, PowerSig delivers competitive or superior accuracy, strong robustness to rough and long inputs, and demonstrated scalability, while using substantially less memory than alternative methods.
4 Conclusion
We introduced PowerSig, method for computing signature kernels of piecewise- linear time series via localized Neumann -oumann-series expansions. Recasting the kernel-defnnig Goursat PDE as Volterra equation yields uniformly convergent, tile-centered power-series expansions of the kemel that propa- gate only boundary data along a directed tile graph and admit efficient adaptive per-tile truncation (Lemma 2.3 Prop. 2.5] Prop.2. 2.8 Prop. A.3). The resulting design achieves linear-in-length memory ((P)), preserves the standard ( nuntime, and supports straightforward parallelism.
Empirically, PowerSig matches or exceeds state-of-the-ar PDE- and DP-based solvers in accuracy, remains stable on highly irregular (low-Hurst) inputs, and scales path lengths previously infea- sible on commodity GPUs. On downstream tasks it delivers competitive or improved predictive performance at substantially lower memory cost.
Future work includes tighter extraction of tile- Mle-boundary coefficients, adaptive scheduling across tiles to leverage hardware concurrency, and refined truncation policies guided by local roughness. Additional directions include extending beyond piecewise-linear interpolation (e.g.. to higher-order segments learned segment maps), integrating certified posteriori error control, and broadening applications in finance, sensing, and long-horizon sequence modeling.
10

---

Acknowledgements
The authors thank Csaba Tóth for helpful discussions benchmarking signature kernel downstream tasks. The authors are also grateful to the program chair and three anonymous reviewers for their helpful and constructive comments and suggestions. A.S. acknowledges funding from the Bavarian State Ministry of Sciences and the Arts in the framework of the bidt Graduate Center for Postdocs.
References
[1] George Adomian. A new approach to nonlinear partial differential equations. Journal of Mathematical Analysis and Applications, 102(2):420-434, 1984.
[2] George Adomian. Solving frontier problems of physics: the decomposition method, volume 60. Springer Science & Business Media, 1994.
[3] Jamshad Ahmad and Mariyam Mushtaq. Exact solution of linear and non-linear Goursat problems. Universal Journal of Computational Mathematics, 3:14-17, 2015.
[4] Fouzia Birem, Aissa Boulmerka, Hafida Laib, and Chahinaz Hennous. Goursat problem in hyperbolic partial differential equations with variable coefficients solved by Taylor collocation method. Iranian Journal of Numerical Analysis and Optimization, 14(2):613-63) 2024.
[5] Thomas Cass, Francesco Piatti, and Jeffrey Pei. Numerical Schemes for Signature Kernels. arXiv preprint arXiv:2502.0847 2025.
[6] Uya Chevyrev and Harald Oberhauser. Signature moments characterize laws of stochastic processes. Journal of Machine Learning Research, 23(176):1-42 2022.
[7] Heinz Werner Engl. Integralgleichungen. Verlager-Verlag... 2013.
[8] Peter K Friz and Nicolas B Victoir. Multidimennsioaa stochastic processes rough paths: theory and applications, volume 120. Cambridge University Press, 2010.
[9] Ben Hambly and Terry Lyons. Uniqueness for the signature of a path of bounded variation and the reduced path group. Annals of Mathematics, pages 109-167, 2010.
[10] Franz Király and Harald Oberhauser. Kernels for sequentially ordered data. Journal of Machine Learning Research, 20(31):1 2019.
[11] Darrick Lee and Harald Oberhauser. The Signature Kernel. arXiv preprint arXiv:2305.04625. 2023.
[12] Terry Lyons. Rough paths, signatures and the modelling of functions on streams. arXiv preprint 1oXiv:1005.4555 2014.
[13] Terry Lyons and Andrew D McLeod. Signature methods in machine learning. arXiv preprint arXiv:2206. 14674, 2022.
[14] Terry J Lyons. Differential equations driven by rough signals. Revista Matemática Iheroameri- cana, 14(2):215-310, 1998.
[15] Sean McKee, Tao Tang, and Teresa Diogo An Euler-type method for two-dimensional Volterra integral equations of the first kind. IMA Journal Numerical Analysis, 20(3):423-440 2000.
[16] Michael Reed and Barry Simon. I: Functional analysis, volume 1. Academic press, 1981.
[17] Cristopher Salvi, Thomas Cass, James Foster, Terry Lyons, and Weixin Yang. The Signature Kernel is the solution of a Goursat PDE. SIAM Journal on Mathematics of Data Science, 3(3): 873-899, 2021.
[18] Amol Sasane. A friendly approach to functional analysis. World Scientific, 2017.
[19] Csaba Tóth, Danilo Dela Cruz, and Harald Oberhauser. A User's Guide to KSig: GPU- 2025.
Accelerated Computation of the Signature Kernel. arXiv preprint
11

---

[20] Csaba Töth, Harald Oberhauser, and Zoltán Szabó. Random fourier signature features. SIAM Journal on Mathematics of Data Science, 7(1):329-354. 2025.
[21] Yuxuan Wang, Haixu Wu, Jiaxiang Dong, Yong Liu, Mingsheng Long, and Jianmin Wang. Deep time series models: A comprehensive survey and benchmark. arXiv preprint arXiv22407.13277 2024.
[22] Abdul-Majid Wazwaz. The decomposition method for approximate solution of the Goursat problem. Applied Mathematics and Computation, 69(2-3):299-311, 1995.
A Mathematical Proofs
A.1 Proof of Lemma
Proof of Lemma 2.3 By definition of r((a), need show that the spectrum of T. is o (Te) (0), i.e., that zero is the only element the spectrum of Te- This is equivalent to establishing that for all À (O), the operator Aid Te is bijection with bounded inverse. Here, id denotes the identity C(O0, 1 To show this, note that it suffices to prove that
Xid Te C(T) -())
a bijection: since (()) is a Banach space, the inverse mapping theorem implies that in this case, the inverse is bounded operator. Thus, is left show that
(Aid - Te) K =
(27)
is uniquely solvable in C(T) for all À E CXO) and all E C(T).
One way to establish this result via the Picard iteration, through proving that (see (15) the Picard iterates Kn given recursively by
(s, F(s,t)+ + 1.. elu, v) v) (u, v) du Ko(s, t) 0,
form Cauchy sequence in (C(1),11-113)
(x,+ (s, - |F (at0) (s,t)
where 3 = (2 1/2
This guarantees the existence of the limit K, = lim Kn, as an element in C(T). Since each K bounded by
K2 (8. [K,()))) Ko(s,()l
k((,,1)
1=0
2 || B(o+))
< Fl ((++))
k=
which is an integrable function on I, it follows by dominated convergence that
(s,)tt) = lim (s, F(s,t) + 1/ lim e(u, K) v) du no0
no0
= F({ + [-f e(u,,)) v) K.,Kn, v) dudu,
so that K, solves the integral equation (27). As noted in [5), uniqueness can be established with another proof by induction. We give precise argument for the convenience of the reader.
For this, suppose that there exists another solution K E C(J) to (27) different from K., Then,
(8,1(8, ((s,t)= = [ elu, v) (K,(u, v) K(u, v)) du de.
(28)
12

---

As suggested in
15 showing
(s, (s, t) 2 || ((++))
sufficient for establishing uniqueness of the solution of (27) For the proof by induction, note that the case = holds trivially since Ko = 0. The induction step n + n + is also straightforward:
[k++1(s,)) 1) - K(s, t) = 1/% o(u, (Kn(u, - K(u, u)) du du
$ llell 2 ||| [// ( el(ute) du dv,
le x 2 || oo g² (E - 1) (e (
(m-1) ||||| eB(ste).
where the first line is an application of (28) Therefore, in C(I), a solution to (27) always exists and unique. Together with employing the inverse mapping theorem, this concludes that the spectral radius of Tg is zero.
A.2 Proof of Proposition 2.2
Proof. We can reformulate as the equivalent integral equation (4) which in turn is equivalent to (id Tpn Po.yy
K ug
for id the identity on C((0, 113) and un the constant one-function u4p = 1 and the integral operator
Tp..v : (C) (1-1==)) - (C([0, 113), l1-1 f f (0.,)) [.S Parulu, v) f(u, v) du
The operator Ton., clearly linear and bounded, and has spectral radius zero, r(Tpp.)) = 0, by
1// = r(T = 0,
Lemma 2.3 Consequently (ef. [16] Theorem VI.6)), we have lim
I Peiv
implying that
bq E (0 3ng E N. !|THH Pax < d", Vn 2 nlg
Consequently, the Neumann series Σn= T is || |l-convergent and, hence,
(id )pex) is invertible with (id T po This, however, implies that K is indeed the only solution of additionally satisfying
T
K = (id T Po. -
T Pa.y
Remark A.1. Note that the n-th Picard iterate is related to the Neumann series via
Kn
ΣT 1=0
Ki.
Pa.y
We remark that this convergence result is independent of any bound on Pz.y and extends a classical result for Volterra equations in one dimension [7] Theorem 4.1]. In particular, Tpe, need not be
a contractive mapping. Repetition of the above argument similarly results in the converegence of
Ene k- -n-n--m " which corresponds to the fixed point iteration for solving (k - Tpm. K f. for any l: * 0 and any f e C(I).
13

---

A.3 Proof of Lemma 2.4
Ain Proof. The (bivariate) power series 11 (s.t) + Ei=o
(1112 s't converges uniformly absolutely
on 00 112 S T1,1 since E
Pia
| I elenal for all (s,t) [0. In particular, 11 (ily
partially differentiable with mixed derivatives (0,,, = PI,1 Σ. p Gi)
Pa.y (s. (s, t), for all interior points (s, of Ti,1 Thus, 11,1 solves the boundary value prob- lem (2) on the tile 11,1, as does K. By the uniqueness of solutions to (2). we conclude that k1.l1.11 = K which establishes 13
A.4 Proof of Proposition 2.5
Proof. Let us note first that any point (s,t) E [0, 1J2 can be decomposed as
(s. t) ($. + (dRa()) T((s))
with () E fin and position indices (k(s),((()) = (a(e-111 1)J +1 (() E)) determined by the location of (s, 1) within the tiling (Th,).. Hence, and by the Volterra identity 4
K (s, = (()))))))) + K(8, maa) K (()(()) m(O) [ S PR(s) K (u, du
(29)
for all (s,t) [0, 1]3. Defining the 'boundary maps' Ikd Tu,1 + (k,, 1....,, 1--1--1b- 1)
(s, t) = K(ok, 1) + K(s, K(),, (s,t) E Te.
(30)
the identity (29) can be expressed as
K(s,t)= (sss). + (M) (s,t) E [].
(31)
Since for each (s, t) Te. it holds that (k(s), = 1k. equation (31 further equivalent to the following (As l)-indexed system of identities in C(T,.)),
(id TK, Th. (k,
1).
(32)
From Lemma 2.3 and basic operator theory (ef. [88 Thm. 2.9]). we know the identities 332 are
invertible for Kg.r and the respective inverse operators can be written as Neumann series in 7 T
Kk,1 (id Tk.1)
Σ Th Tlnn
(33)
where the above series converges wrt. Illoo: the sup-norm on C(Ti.)).
For the recursive structure of the Ka.r-identities (33, note that, for any fixed (s, E Tk. we have
(ax, 1) E T n Ta, and (s, π) E Tk. n k. where by definition To {0) [T. 7+ i] and
Tk (0I, 0411 {0} and To,0 == {00, 0)). Consequently, the boundary map (30) evaluates to
(s. 11-1, t) K (s, m) FR-1,1-1 Tt),
(s, E Te
(34)
with Ko,1 = K To. Kk.0 = K| Kl K0.0 = 1. Combining 34 and (33
proves 15
Lemma A.2. Let y0 1]2 (s, + (s oy and ot [0, 112 (s,t) + (t - T)', for any given 1 E No. Then for each No and all p. v E {1, (1.....
1}. we have that
Piuy
T (v (T.. (8, = (1+ 1)
(s kp) + Tz)", (s, E T Taws and
u.v T
Piu.
Tm Tmp((p!)) Tp.. (s,
(s -op) (t Tx) ,ltt (s,t) Tp
(1 + 1)
14

---

Proof of Lemma A.2 for n E N we get
This follows immediately by induction. Indeed: the case n = 0 is clear, and
( T Tun (s,t)= = Pa.v
Thi T
(1) ( Nu, v) dudu
op
Pp
(n (n-1)) (u on) (w)(tn-((w =u)m² du de
Pa.
(1 + 1)m-1(nn
ou
PM
(0- (w- (i-To)--dhe de / (u- 1p) du Jop
(1 + 1)m-1111 (n 1)!
Pp
(1- Ty) (s 0 ptn
for each (s, E T Tue:
(1+ 1)m- 1)!
1
as claimed. The proof for TH (p! !1) I is entirely analogous.
(To.)
A.5 Proof of Proposition 2.7
Proof. We proceed by induction on the tile position (k, For this, note first that, for all o, 2 0,
L = (ō.g - o o20x0)
(1
k=== Ei 1in-a6600 600 )8k jōk-1.) 1,520
1.120
A10
1,20 = (lir(1))) and
(35)
j+1
Rs (1.1-1
(6i5 - 7'00))
1 (Enmo
(2z - yos 1.520
i.20
bi. 10--1-111:0 1 6os
1.50 = (rag(T))is20-
(36)
Then for the base case (ke, 1) (1, we have that (C
C...
C0,00 (a,t)
1 and hence)
C... (ōur 60z ).520 and further, by (2I). that
Lo Ro = (ei;(0)) (En= (doi 8an)r,;(0) 1,220 = (l:)(0))s520) (doiTo) (0)) (0))
= (En= in (0)ros(() (0) don).,20 = (Coo(O)ro)0)), i,j20 = (5 61.5) 1,520
Hence [cl (a.r)
[pI.] Lo Rol (a.r) = -1,,1,,, whence
20
(k.1)=(1,1) holds by Lemma 2.4
(and via induction n) if, for any fixed n E N22.
p (nly
den dns)
1.20 = -1,1
Pi [jn=1 dnn can-16n-11 18n
Ro.
(37)
i,20
Denoting aij
[[n [[n-1]]j= i, - 8n- ij for each i,j 2 O, the right-hand side of 37 reads
Pi.
P1,1 Lo (E - Unn T us (0)) 1220 P1,1 (E (o anj (0) i.20
(38)
Now by 35 and 36
we obtain for any fixed i, U, H 2 0 that
L (0) in (0) Olyp
8, 1, 1 brut1 @-w=1 16-11-11
*+1
1-1
This allows us to evaluate 38 to
Pl,1 (Zve (0)
(0) un ug (0)) 1.500 PI.1 (ai- 122000000000000 P1,1 (
(Im-ij)pps din
P
dindns iss20'
proving (37p as desired.
2-112 With the base case (k; = (1, 1) thus established, take now any (k, E {l...., f- 1]2 with
k ++Z / 2, and suppose that (20) holds for (k: (,,) and (k: 1) and (k 1,1 1)) (induction hypothesis). Then Kg.i
((a.r)
for each (k. E {(k: 1, (k: 1), (k. 1,1-1 1)}, whence (a.7) Te.r
and by Proposition
2.5 we have that
Σ n=0
(* + (n(-), ((())))((m())))))))))))()))))( + (nr(ax), -*+ (39)
KR.I
15

---

uniformly on Ta. Since Ty = id )(T,))) the 0th summand in (39 reads
(((o)) + (n))), C ((()) + (n(ox), C m(()) [
(a:))
for the initial matrix of from 22).
Consequently, the claim (20p follows if, for any fixed n E N,
chid n+1
-Tk, (1 ((rr) loa.).
(40)
(ast)
Th.1
Abbreviating un+1 = Thi [Chit] (omr)
and Ofe = (cty)s note that by definition,
1 [ s'E d = Pk.l Σ
un+1(s,t)= Pk.t Σ
(i + 1)(j +1 (5 - )f 7/
Cjj
Cij
71 nov
i220
1,20
(41)
for each (s,t)= E I (see
(14 and 19 Abbreviating (Ciz) = ((++1(1++++) (++1)(+1) s.520 and using that
- ( (ft -//1+
s 1+11 o 1 - 7/+1 gitl ok +
for all i.j 2 0, we can immediately rewrite 44 as uk+11 = Ö(ar)) nv.
for the coefficient matrix
Č == (C 11-11: dos M/-11 510 + 9:: 520005 1,200' ),200
where 7gl: -i= ب and 7:15 Eizo k
ftl (i.j 2 1) and 9-11: -:111 Ć-1,J
- 0 (i.j 2 0), and í4:: Ei.20 i+ 7f Hence, the desired identity (40) follows if
Č = Pho Chud R
(42)
For a proof of 42 note simply that, by definition of Lo and R, (ef. 35 and (36
Lop Chid Rr,
(oxk) 1.520
(Ennzoldrwtt of u+t dea) [u+l]lo+1) (6ms-1-10) - +i +1-1160)) 1,20
Cur
(ci-1--1 E 1i 1 dos -p- ooth 6i0 + Emozn u+1 7 dindos 1,520
(Q-15-1 8/11 doz h---1 8.0 + t..: 6zn0o;) 4,20'
This implies
42 and hence concludes the overall proof of the proposition.
A.6 Proof of Proposition 2.8
Proof. We combine Lemma A.2 with a similar reasoning as for Proposition 2.7 Proceeding by induction on == k + I, we will show that, uniformly in (8. () E Tx.
(s, TI) ( (R,0) (s ox)" and (k,,(o =
(kh,)) 77
(43)
for each E (1,...,-- 1-1 (which, particular, implies that an
(k:0) 1lo Ti) = g(kdd)) and
use this to then prove (24). Starting with the base case 2, note that (43) (1.1) is immediate
from the boundary conditions in 2 and the definition of o
(1.1) and (( Thus, for (k,1) (1, 1),
(1)(
+
(h.)) TT (v(f))
(44)
3
16

---

by combination of 43
with Proposition 2.5 which holds uniformly on Tk.1 and for the zero-
starting sequence
$20 = ((1 - 6, (k,d))
$20" Denoting Un = Ei a. T,TT,
Ej=o B, B, (4pt)
for each n E No, we obtain that, pointwise for each (s.t) E Ta.s.
| (5-oz)) ow) d8 f(i (((--))a n)
(1,1( = (Tr,(un))) PR.
Σ (k.l)
1.==D
Σ (s oR)"
i+1
() + +) (t- (t-n)) nTttt = (n1s oox): ( I)n(t Ti))
+==
(45)
where from the second equality onwards we assumed that = (n( - ow), Oh.- T() for some (c()dtm) l (N8) (induction hypothesis). Since indeed
Ch (c
20
Ug Σ
(kid) 0i
(s ok)² + 1.. - TI) = (vo( vo(--ok)) ok), C/ -1)
(46)
for Akll = (a,
(1.10 (/.1 810 + }3 g/ oos 1k-1,1-1 71)6u: d0z )0y) combining 46 45
proves
Kk.1 (n( ox), ox), Cl --))
for Ch. Σ
Ik S Cp ck.
(47)
n=0
by induction on n E Ng- (Note that at this point, 47 is established for (k, = (1, only.)
Since (5")
516. and (I")ss ([+j)!
it (i-nj! for each i,j 2 0 and E Ny by definition of S. T.
n and j 2 n,
n!y!
S S C.k.
a (1 (i-
(k.))
1.j20
in!
2 n and j = n,
0,
else.
Consequently, the matrix entries of Stiil =: (c
(k.1)
1.,20 each read
min(sal) PEd g (i-min(t.)))
2 min(i.5)
mints min(i Ped a (i-min(k,))! mim)
(Ar. 2 Be, W)
(48)
i min(i,j)
implying G C..., desired. But then, in particular recalling 47 and the definitions 26
Kk,141 T(11) (n(s σk), Ot c - 71) -n1)
(n(s on), (a (+.+++1)) )+20) (
(1,11) -O and analogously
(0H1)(00111)) ((3(6+1.)) m20, nll (((0-m)))ees
g(++1.1) (1-
TT)", (1-T)"
uniformly in (8,t) E Tk. This-together with the fact that Kk, (8, 1)) = 11 Di or (s ox)" (k.1)
(k.1)
for (alk::) 2r20 (80,)) and all [o ok+]] for each k 1,...,f--1 1) and K1 (0,1 1,, t) = (= Eee ((1.)) (t T1"" for (e 2, = (8u,))22 and all t E [0, 11, for each | 1.
1) implies that 43 holds also for 7 3. fact, the above argumentation-aa stated-proves that
stated
if (43) holds for some fixed k, E [1....,f -...,,-11 then both (47) and (48) hold for this (k, and
also (provided max(k, 1) S - 2) that (43) holds for (k: + 1,7) and (k., + 1). This proves the proposition.
17

---

A.7 Gram Matrix Approximation Error
The remarks on local truncation error in Section C.3 can be directly extended to yield an explicit prior bound on the approximation error for the Gram matrix of given family of time series: the doublc-sum structure of the local approximation error (59 allows for direct control through modified Bessel tails, which can then be scaled from component-wis to matrix-level bounds.
Note that an estimate very similar to the one underlying the following proof was first presented in S
Proposition A.3. Let G = (Gis) = (Kat0, a0 (1, 1)) b.j= be the signature-kernel Gram matrix for a family X = (xx(1 (a) Ii=...... e 1 1, ,,...,m) of time series in R. Let further where 2! R is defined as (24) but for the matrices
GN = (N(N) (1 wil1 (1,1) ( i
CN = A B WN
(k, (1....,,-11]
min((s)) with A!7 = (PEN
OSIISN and WN == (wis)oSissN and B = (ViJ) (D.NDIN))
J for
(R.D)N
= a ((.D)N
((KD)I and 6 M/r
for each (i,r) E E
1stt
where B111IIN = or (2,1)1A = 5p. for each r E No and, recursively (recalling 26 for notation),
(a.K. N -=0 = (k,(1-1),
(n(((((e-1)) (A,(L-D.,(I/(E--))
and
==
D. m(1((e-1)) (1-(())
for each k, E (1.....-11. Then we have the (truncation induced) approximation error estimate
1 CN
c 11.
(49)
(N + 1)!2
with || the Frobenius norm, || == maxisuujsmisssstt xxx Aاz(امي and with
26-2
N+1
II 1o
2 Xl
and SN = [
2/-2
E-1
N+2
8-1
Proof. We proceed along the aforementioned lines, starting with the trivial norm inequality
|G- S m
max Sm
Eij
[N]
for Eg = Knti, min
[N] min (1, 1-------11 1: 1000 (0 (1, (50)
Now, for any fixed (i,j)= E {l...., m]x2 x2 the error Sij reads, see (the last display of) Section A.6.
[N]
IN Ejj (1,
IN ( Σ El (1-
E6 oe-1""
(51)
for some coefficient sequences ge'l' = (Ef.. )r20 and Ee, = (EC )r20 (where we suppressed the dependence on (i,j) and N and denoted E' = 1 to ease notation). By linearity of the Goursat PDE (2) and its uniqueness-of-solutton
we find by comparing coefficients (see Section A.6) that:
ct,'"' L. -ath" and ge'L =gee -30h"" 30'.'' for the formerly defined coefficients a. (af )ro
(k.))
and any =a (K.)JN 0r20 and ghd (8)
(g(k))) Drzu and gli.d (af (k(k.))) Drz0. and that (as we recall)
(a gk.d) 24 (akk- gk-111 and (a (44))
for the (bounded) linear 'ADM-type' Goursat-solution operators 24. Ax. (N)) x2 e (No) x2 defined (recursively) in Proposition and and (via projection & (1%) - R2 N Proposition A.3 respec-
18

---

tively. In this conceptualization, the coefficients 15 of the approximation error 51 read:
(5 (gle = (o -0" - Ue, (a 18
(M(N.o)) (a -- + (e)e) -0 Aee
(a (a/ -1 1,00-11 1,
= -T((NoO) o (a 6'''-1 -1 + (T. o2lgr (a - at/ -1 ge -1.' - -1/ 36--1107
(T(N. (a (. 10000-000 -10) + (To.No (1 -5t-
- (Zl (a 11' 3/- - + (TONN)
6 -1 ,² 11.00)
(52)
(TO.NO ( -1 --t -1,0 1,00)NN
for the projection Noo : (No)) x2 a (af" a 120 + (of Lo. (r), (2) No.
(1)
ASN E (M and its defect T(N,=)) = (16) - To..N| with a = 2o)) (a) for each a E f (No). Note that the third of the above identities holds by the linearity of 2,, and since the
operators N] Ila, and I,,, coincide on the subspace (T (No) x2).
Denoting the summands in (22 by 5(),, 5(2), and 5( resp., direct (but tedious) computations show:
max 1
(1) 11 S (2(2" - (!2(e')
(53)
2 2//22' 1)1lvl
(26 (
F!2 (1"(("))1EF)
v=1,2
max +
(2) 1 S To
8-1
2((')r122
E=E Me E 1220 ((28' -1) l)F1
max
(54)
max
(3) rl कe N Io
1)1ssal 2/ 1-1
(28" 2(e')
max
(()((())
(55)
g ,80,0 ((28' 11)xx,
1 S (((2" -2) (122" 1)
(n=
(56)
r!2(e))r
for each E No. with Iu : x + Io(a)) the modified Bessel function of the first kind of order ( ماz(111, (,2(")) and for the function pp() =
zero (cf. Lemma 2.4p and [F5] == maxigk.Se-1 11 1,,
II To / Fig (()) and NN = (PyCa + N N 1
see Sections A.5 and A.6 and ef. also
5 Props. 3.2, 3.3, 3.4). Applying the estimates (3).
555 (56) (on the size of the additive
components c(0) in 52 back to (51) via the triangle inequality and the auxiliary estimates
(20 (ea))
1 (
Σ 21
S 00.5 10
((2 11)ee,N .N
1,1
E-1
(g)r!
(20 120") D"
26'11
(28') N+2
(22
S x(22 1)N
Σ n N+1
S p(22 1))eeN
(r Nu)
(6')r!
N +
(e')r!?
n=0
(57)
(E" (126-1)]Bs) -1]11s Tr and NerN = 126' + then yields the error bound
with 05.5 == max g C gf ,20
[N]
El + Err
29' ((')
Kllll' + 26'26 ((0)==
(58)
N=
N
= Σ
(((2) + Sp
(3)
1,1 (e') - +
b. (C)
==
+1
Np.N (2 N+2 (
00
(28 k D"
S 522 -1)
Σ
N +
(0)2rr!12
(g))2y!
2N+2 [zy) N+1
2N+1 |E N-1
1k1/ N+1
S ((28))
(((2))(N
(P) + 1)N +2)
(c))+1 (N 1)!²
(N + 1).²
19

---

Note that the second and third inequality in (57) follow from (54) and (55) via a straightforward in-
duction, see also [5] Prop. 3.3]. The desired inequality (49) now follows immediately by combination
of 8 and (50). using also that the function [0, oo) Io /ypu
is monotone.
B
Benchmarking
B.1 Downstream Experiments: Additional Figures
Complete implementation details and hyperparameter grids for the experiments underlying the following figures (see Section 3) are available the PowerSig GitHub repository at
https:// http://......eeebbbtt///rrrrrree
Training
Two Day Rolling Average
Towerds
17500
12300
10000
2000
102
150
299
Day Rolling Average
Actual Peweriic (4415)2.815)
8500
1000
15
25
30
Figure 4: Bitcoin price regression (two-day rolling average). Top: training fit; bottom: test fit. PowerSig attains 2.81% MAPE (under the default linear static kernel) compared 3.23% MAPE for the RBF-assisted KSig-PDE, while using only fraction of the device memory.
20

---

Eigenworm Classification: Accuracy vs Window Size
KSigPDE KSig RYSF-TRP
PowerSig
0.0
0.6
Accuracy.
0.2
0.0
Window Size (log scale)
Figure 5: UEA Eigenworms classification across window lengths L. Test accuracy versus input window length L for PowerSig. KSig-PDE and the RFF baseline RFSF-TRP. PowerSig remains competitive and scales to L = 1048 with 61.1% accuracy; KSig-PDE is competitive to L 128 before running out of memory (OOM). RFSF-TRE TRP attains slightly higher peak of 62.5% but OOMs for larger L. consistent with the storage advantages in Fig. 3 This shows that extending the input window, enabled here at scale by PowerSig, can narrow pertormance gaps often ascribed to inductive bias while preserving feasibility.
128
Timesteps
Timestaps
87S-T33 Powerflg
FFST-19
367
boar Abouk ta Nercer Cage Entrr (MAPE) No
Fower
10'
Moon Squa Erorr (M66)
0012
10
1077
1000
1500
1750
2202
3250
3000
3750
1000
1258
1585
1750
3000
2250
2000
2750
Number of Timestegs
Number of Timestaps
Figure On long-horizon periodic signals, SVM-regression error versus input-window length for synthetic pear-periodic time series with adjustable period. For representative instances, error decreases monotonically as windows span multiple periods. PowerSig sustains this favorable trend at window lengths beyond the reach of conventional or low-rank signature-kerdd baselines, while maintaining peak memory that grows linearly with window length.
21

---

Dimension
-9.5
-1.0
3300
3600
1960
3850
2000
Sample Dimension (Trrget +1
1.0
0.5
-1.2
1300 1460 1500
3600
1760
1800 1900 2000
Time Stup
Figure 7: Representative near-periodic synthetic time series with adjustable period length (indus- trial/sensing proxy) used in the long-horizon experiments above. This instance illustrates the quasi- periodic structure across multiple cycles on which SVM-regressio is evaluated; as input windows for SVM-regrsssors span several periods, their error decreases monotonically, with PowerSig sustaining this trend longer windows while preserving linear peak- memory growth (Figure
B.2 Comparison Against polysigkernel
We benchmarked and contrasted our method (PowerSig) against the concurrent work polysigkerne from Cass et al. [5]. Both our approach and theirs were developed independently and released within two weeks of each other, with our preprint and theirs sharing the same arXiv month stamp.
Using the public JAX implementation of polysigkernel, we ran identical experiments on com- putation time and memory usage over two standard benchmarks: (a) fixed-length (= 512) paths of increasing dimension, ranging from 2 to 4096, and (3) 2D Brownian motion paths of length 512. The results, shown in Figures SaaressSane and respectively, show that both solvers achieve essentially identical accuracy throughout while PowerSig performs better runtime and memory usage. In benchmark (3), PowerSig was on average 3.1% faster and required 23% less memory than polysigkernel on oscillatory time series of length 512 (averaged over 10 iid samples). In benchmark (a)) PowerSig exhibited a comparable computational advantage, being 8.5% faster on average at dimension 512, with similar memory usage.
Runime Dimeniio 4298)
SJU Memory Usage Dimenson 4298)
Powet5p nadckplanel
TowerSig MySkplaam
107
Oiseniion
Chrandon
Figure 8: Scaling wrt. dimension (benchmark (a)): Runtime and peak memory fixed length 512 with dimension ranging from d = 2 to d 4096. PowerSig is 8.5% faster on average at dimension 512, with similar memory usage.
22

---

Accurscy Comoarison: Poweriip vs PotySig
Runtime Compartson: PowerSig PolySig
0077
103
M406 (woatoo (la)
10
Figure 9: Scaling wrt. time series length (benchmark (3)): Runtime and peak memory on 2D Brownian paths with length ranging from f = 2 to L 512 for PowerSig and polysigkernel. Accuracy is essentially identical, while PowerSig is 3.1% faster on average and uses 23% less memory at length 512 (mean over 10 iid samples).
Let us clarify why these empirical differences occur and how our underlying methods differ in concept despite propagating the same local power-series map: polysigkernel starts from explicit Rinntnn-function formula to rewrite the Goursat solution as sum of integrated modified Bessel functions [ Thm. 3.1 (Polyanin)] which, after expanding each in its well-known power series, allows the kernel on every rectangle to be reduced (using -eeggrationbb--aarr and change of variables) to univariate centred power series in s and for which the resulting coefficient arrays can then be generated by the closed-form recurrence ACila..ee x(T.),.1] in [5] Props. 3.1, A. A.2]. Truncating these expansions up some fixed degree yields recursive polynomial approximation scheme mollssigkerner for which the authors provide local and global truncation error bounds [5] Prop. 3.4: Thm. 3.2]. power-performan
scheme for polynomial boundary data is also discussed in [5] Section 3.3].) Our approach, by contrast, approaches the same PDE from a functional-analytiic view: Using the Volterra formulation of the Goursat problem, we express the kemnel propagator through an integral operator of spectral radius zero (Lemma 2.3) to obtain concise, three-step Neumann recursion (Proposition 2.5) whose iterates build the same local power series ceffficienttby-cofffcciieee in way that allows for their immediate coefficient extraction on the fly (Proposition (28). The resulting iteration adapts per tile and stops once the next term below machine precision (typically 5-8 iterations; global cut-offf N required), leading to one-strip memory profile that below the quadratic footprint of the reference polysigkernel implementation.
ceebyshev-interolattee
In summary, polysigkerne obtains each tile's coefficients from closed-form Bessel (or Cheby- shev) recurrence evaluated up to user-chosen degree N, whereas PowerSig builds the same coefficients adaptively via a Neumann iteration that stops once the next term drops below machine precision. Both implement the same local propagation map and are analytically interchangeable per tile, with their different constructions yielding distinct practical profiles with modest but consistent performance edge for PowerSig our matched head-to-head benchmarks (Figures and g)
C Summary of PowerSig: Methodology, Implementation, and Limitations
C.1
Algorithm
The core idea behind our proposed signature kernel approximation method, PowerSig, is illustrated in Figure of the main text. summarizing deseription, in pseudocode, of an efficient Python implementation of this method-ale used for our benchmarks is provided in Figure T0 below.
C.2 On Benchmarking
We benchmarked algorithms systematically terms of runtime, memory usage, and accuracy. To ensure consistency and reproducibility. we developed standardized benchmarking framework
23

---

Algorithm: PowerSig
procedure PRoCESSDACCOAASS(T(...ooo co. r rows-1] order)
2: Initialize ( 00] < 1,0,...,0
Initialize (0)) 0 <
0
4: for d + to rows + cols 3 do
(cols 1)) end_row 4 min(d, rows 1)
start_row max(0,
6:
7:
en-eworow- start_tattoo +
8:
Initialize s gl L), ( [0.. L] zero vectors
9:
for 4- to L 1 do
10:
i E sttattttroooo row + k
11:
j d-d-i
12:
if rows and + 1 cols then
13:
Compute p + (X(j + 1] X j.
Yl 1] - Y(i)
14:
Form Toeplitz matrix U from (d) [k] (first column) and s(4) [k][1 :] (first row)
15:
Form Rıj 5 min(t..) for i, o..... order
16:
Veol f [a ,......."" under??
a cols
17:
Vrow [30
gorder]
B
18:
Snew Viww (U R)
19:
nnww
(U o R) Veol
20:
if j < cols 1 then
> Before right edge
21:
if k > 0 then
22:
s(d+1) + 1] f- Snew
23:
else
s s(d+1)(0) [1.0...... 0]
Initial condition
25:
end if
[d+1]((k] tnow
27:
else
At or after right edge
28:
s [d+1] 1) [k] f- Snew
29:
if k 0 then
30:
(d+1) 1] +- fnew
else
32:
t [d+1)[kk] [k] + [1,0...., 0]
Initial condition at far edge
end if
34:
end end if end for
36:
37:
end for
38:
end procedure
Figure 10: Our implementation PowerSig, g, shown here, is based on the algorithmic approach that we developed in Section 2 and summarized in Figure] ]
24

---

comprising a base benchmarking class, a custom context manager, and a proxy wrapper of the CuPy allocator. While this setup directly supported Py PyTorch and CuPy, special considerations were necessary for JAX, as its XIA backend does not natively allow resetting peak memory usage between runs. We circumvented this limitation by isolating each JAX benchmarking run in separate Python subprocess, ensuring accurate GPU memory measurements and preventing JAX's s default behavior of pre-allocating 75% of GPU memory. Additionally, to prioritize numerical accuracy required by our experiments, configured JAX explicitly to use 64-bit (f64) floating-poin precision, overriding its default preference for computational speed.
C.3 Error Analysis and Limitations
We identify two primary limitations of our proposed method for computing signature kernels: numerical stability (truncation errors) for very large Prry (see also Proposition A.3 for details), and computational overhead due to JAX recompilation.
Our computation of the signature kernel effectively involves linear operations to evaluate polynomial expansions, the latter given by the power series (24) for which by construction:
Truncation Error
min(t.5)
(s (t- -int-Ti TiY
[N| Nk.
(s,1)
Σ
(s (t Ci
- mI
PR.
((s, E Tw.)),
1,j=0
m==ltj-minn pV miim-m(()..))
for any truncation order N E N. Evaluating this at the tile boundary point 74141,7+1) then yields:
DC
[N] K
I min((t.f)
(Ao)
(59)
TO:
i,==+41
for Ao = 0-+1- =0 o% and AT == T+- As Ao, AT 0 with increasing time series length, terms involving (Do)' and (Ar) will rapidly fall below machine precision a speed counterbalanced by the magnitude of the |px.1l values; large values of [Pt.1| significantly amplify the truncation errors (59). necessitating higher polynomial orders N (which again are bounded by numerical precision).
However, truncation errors are typically localized, decaying rapidly across tiles, which helps to limit the overall impact of few tiles with higher error the accuracy of the final solution. To further mitigate truncation issues, our implementation configures JAX to use 64-bit floating-poin precision (f64) and includes a mechanism for estimating the minimal required truncation order N meet specified truncation error tolerance. Roughly speaking, we evaluate the intermediate ADM coefficient matrix for the maximum |Pk.1| in the grid and then compute the sum of all entries in the last column, scaled by A and in the last row, scaled by A This provides bound on the truncation error and allows for a straightforward search over N = 8, 64 to identify a suitable truncation order for given set of coneet-cnstituting time series.
more detailed discussion of approximation errors, including a rigorous analysis of the local truncation errors (59 is provided in Section A.7
JAX Recompilation JAX's compilation mechanism can incur significant overhead when pro- cessing batches of time series with varying lengths ("jagged" datasets). One potential mitigation strategy- using re-interpolating to pad each series to a uniform length-in not ideal, as it may distort the underlying temporal spacing and thus affect signature kernel computations. However, most datasets are trimmed and chunked to the same length so this is typically a minor limitation as long cardinality of the set of input shapes is not approximately the same as the number of time series being compared.
25

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification:
The abstract and introduction include the claims made in the paper.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
It is fine to include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: the appendix (see supplementary material), we address truncation error and JAX recompilation potential limitations in applying our approach.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
the paper. The authors should use their best
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
26

---

The theoretical results of the paper are presented as fully proven lemmas and propositions, with all assumptions explicitly stated as premises of those results. All proofs can be found in Appendix A of the supplementary material.
Justification:
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proof's can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification: Accurately benchmarking signature kernel implementations across diverse frameworks (e.g., NumPy, CuPy, JAX. Torch) terms of runtime, GPU, and CPU memory usage is notably challenging. To address this, our open-source library provides a context- manager-basee based benchmarking framework, featuring abstract base classes that facilitate reproducibility and transparent benchmarking. Additionally, we supply all utilized seeds, runtime parameters, and data assets to ensure reliable evaluation and reproducibility.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language model), releasing of a model checkpoint, other means that appropriate to the research performed.
also be provided via detailed
While NeurlPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make clear how to reproduce that algorithm.
(b) the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully.
(c)
If the contribution is new model (c.g., a large language model), then there should either be way to access this model for reproducing the results or a way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility.
27

---

In the case of closed-source models, it may be that access to the model is limited in some way (e.g., registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, described in supplemental material?
Answer: [Yes]
Justification: See the answer Question 4. All source code and parameters supplied in the supplementary zip folder "powersig-code.zip."
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips.cm cc/ sublic/guidess oowpss//oooesubbisssoooooooooooo
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
//nips. http/uullcggguucee https//oooSSumiissseeee asiionooicyy for more details.
The authors should provide instructions data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: Yes, see the answers to questions 4 and 5. We also include hardware con- figuration and JAX optimization settings a file in code-documenting zip-folder "powersig-code.zip"
Guidelines:
The answer NA means that the paper does include experiments.
The experimental setting should be presented the core the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [NA]
Justification: We benchmarked 10 rounds per algorithm with warm-up round. We tried adding error bars around the mean but they were so tight that they didn't 't render well on the graphs so we removed them as they didn't add much value in contrasting the relative performance of the algorithms.
28

---

Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by enror bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall
with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
should be clear whether the eπror bar is the standard deviation or the standard error of the mean.
It is OK to report I-sigma error bars, but should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not show in tables figures symmetric error bars that would yield results that are out of range (c.g. negative error rates).
error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
Question:
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Answer: [Yes]
Justification:
Yes and we used workstations with high consumer grade GPUs.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't 't make it into the paper).
9. Code of ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics http:p///eerrrpp..c////////eeeie
Answer: [Yes]
Justification:
We have reviewed the NeurIPS Code of Ethics and can confirm our compliance.
Guidelines:
The answer NA means that the authors have not reviewed the NeurlPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g.. there special consid- eration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive socictal impacts and negative societal impacts of the work performed?
Answer: [NA]

---

Justification:
This paper presents foundational research and is not tied to specific applications, let alone deployments, that we are aware of as having any notable societal impact.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No. they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation. generating fake profiles, surveillance), fairness considerations (c.g.. deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (c.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: To the best of our knowledge, the release of our algorithm does not pose any immediate high risk of misuse.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-us, should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification:
Creators or original owners of assets used in the paper are properly credited.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
30

---

The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g.. website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperssi thcode com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that re-packaged, both the original license and the license of the derived asset (if it has changed) should provided.
If this information is not available online, the authors are encouraged to reach out to the asset's creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
The only new asset we introduce is the library for computing signature kernels and the benchmarks themselves. itself only relies on other open-source libraries for implementation and benchmarking.
Justification:
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification:
The paper does not involve crowdsoureing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then much detail as possible should be included the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: The paper does not involve crowdsoureing research with human subjects.
31

---

Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with
human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is important, original, or non-standard component of the core methods this research? Note that the LLM used only for writing, editing, or formatting purposes and does not impact the core methodology. scientific rigorousness, or originality of the research, declaration is not required.
Answer: [NA]
Justification: The core method development in this research does not involve LLMs any important, original, or non-standard components.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, or non standaadd components.
Please refer to our LLM policy for what should or should not be described.
32


Reviews:
[{"review_id": "SjhLtfUYAw", "reviewer": "Reviewer_VWgU", "summary": "The authors present a scalable approach for computing the signature kernel that leverages adaptively truncated local Neumann series expansions of the Goursat PDE (a hyperbolic boundary value problem) associated with the kernel. The authors show that the proposed approach is applicable to time-series with a high degree of roughness/volatility, enables massive parallelism, and achieves orders of magnitude reduction in memory requirements compared to dynamic programming and PDE based approaches.", "questions": "Here are some suggestions for additional numerical studies that I believe will strengthen the submission:\n\nTo make a compelling case for PowerSig, I suggest that the authors provide numerical studies on at least a few learning tasks; e.g., Multivariate UEA Datasets, https://github.com/tgcsaba/KSig  and compare against alternative approaches, Given the efficiency of PowerSig, such studies should be feasible within a reasonable computational budget. This will also provide some insights into whether the benefits of (almost) exact signature kernels translate into superior predictive performance in practice.\n\nNumerical studies on higher-dimensional paths (say d=10, 50) to confirm that the runtime grows linearly in d and accuracy remains stable.\n\nNumerical studies showing how P should scale as a function of roughness and sequence length to maintain a target accuracy.\n\nMinor suggestions for clarity:\n\nRename the baseline explicitly in both caption and axis label, e.g. ``MAPE relative to KSig truncated-signature kernel (level 21).''\n\nKeep the baseline constant across panels or state the new baseline explicitly in the right-hand axis label.", "limitations": "yes", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "The authors demonstrate the ability of their approach to handle time-series with over $10^6$ time-stamps on a single GPU with 24 GB RAM. This work therefore presents an impressive advance over existing methods for near-exact kernel computation over long paths.\n\nThe paper presents a linear sequence of lemmas and propositions that rigorously justify, step-by-step, the key ingredients of their computational strategy. This principled, structured presentation enhances clarity. The authors show that their approach ensures convergence when the truncation order goes to infinity. \n\nThe authors present numerical studies for a two-dimensional Brownian motion test case and compare against KSig PDE. Results are also presented for two-dimensional fractional Brownian motion paths of fixed length across decreasing Hurst indices. The results serve as an empirical demonstration of the memory reductions and runtime scalability for long sequence lengths.\n\nFrom a numerical analysis viewpoint, it would have been more insightful if the authors provided theoretical results establishing bounds on the Gram matrix approximation error as a function of the truncation order. This would also provide useful insights to practioneers. \n\nThe paper does not compare against Random Fourier feature (RFF) approximations or low-rank approximations, stating their tendency to degrade in accuracy for larger time-series length or time series of higher dimension. While this is correct, in practice, accurate kernel computation is necessary but not sufficient for good generalization. Approximations such as RFF introduce their own inductive biases that can be beneficial for some datasets. \n\nAll the numerical studies presented in the paper fix the truncation order P=7. The results show error with respect to a reference KSig truncated signature kernel. For a machine learning venue, additional numerical studies on learning tasks would significantly strengthen the submission.", "was_revised": true, "final_justification": "The authors have thoughtfully addressed all the points raised in my review. They have shared results from additional numerical studies  they conducted to address my concerns. This is an excellent paper and I am happy to raise my score from 4 to 5.", "rebuttal": "Thank you very much for your time and your very valuable comments and insightful suggestions. In what follows, we address each of the points raised in the review separately, hopefully answering and clarifying any open questions or concerns. Kindly let us know if we should elaborate any further on any of the points. \n\n1. From a numerical analysis viewpoint, it would have been more insightful if the authors provided theoretical results establishing bounds on the Gram matrix approximation error as a function of the truncation order. This would also provide useful insights to practioneers.\n    \nReply: Thank you for suggesting this. We agree that an explicit bound on the Gram matrix approximation error would be valuable, particularly for practitioners. To address this, we will add a corresponding result at the end of Section 2. This follows directly by refining our existing remarks on truncation error in Appendix B.3 into a fully specified inequality, which is straightforward: the double-sum structure of the displayed approximation error (as displayed just before (49)) admits direct control via modified Bessel tails, which can then be straightforwardly scaled from component-wise to matrix-level bounds.\n\n2. The paper does not compare against Random Fourier feature (RFF) approximations or low-rank approximations, stating their tendency to degrade in accuracy for larger time-series length or time series of higher dimension. While this is correct, in practice, accurate kernel computation is necessary but not sufficient for good generalization. Approximations such as RFF introduce their own inductive biases that can be beneficial for some datasets.\n    \nReply: Thank you for pointing this out. We are currently implementing a comparison against low-rank kernel approximation schemes by benchmarking $\\texttt{PowerSig}$ against the recent specialized low-rank and RFF-based methods introduced in [Tóth, Oberhauser, and Szabó: “Random Fourier Signature Features,” *SIAM J. Math. Data Sci.*, Vol. 7(1)]. We expect to have results within the coming days and plan to include them in the camera-ready version.\n\nOn a preliminary note, to evaluate performance in a downstream learning task and to compare with RFF-based methods, we obtained preliminary test results on $\\texttt{PowerSig}$, RFSF-TRP, and traditional linear/RBF kernel SVMs on the standard Eigenworms time series classification dataset for increasing time-series lengths, from within the UEA dataset that you kindly referenced. Due to the extensive preprocessing involved in such tasks, it's difficult to fully assess the potential of signature kernels on long time series out of the box at this current moment. However, to highlight practical applicability, we have obtained preliminary benchmarks on accuracy, precision, and recall across multiple truncation lengths of the Eigenworms dataset ($L = 50, 100, 200, 300,$ etc.). Over this range, we observed the accuracy of $\\texttt{KSig}$ and RFSF-TRP drop from 48% to 43.5%, while the accuracy of both $\\texttt{KSigPDE}$ and $\\texttt{PowerSig}$ increased from 32.82% to 40.46%. We also observed that, contrary to $\\texttt{PowerSig}$, RFSF-TRP begins to run out of memory at length 1000 under default settings even on an H100 GPU with 80 GB of RAM.\n\n3. All the numerical studies presented in the paper fix the truncation order P=7. The results show error with respect to a reference KSig truncated signature kernel. For a machine learning venue, additional numerical studies on learning tasks would significantly strengthen the submission.\n\nReply: Thank you for this suggestion, and for the reference dataset provided below. We fully agree that downstream learning tasks would be an important benchmark from a machine learning perspective. Unfortunately, practical constraints have so far limited our ability to include such experiments at a convincing scale. Specifically, we encountered severe resource limitations when attempting to reproduce kernelized learning tasks on the full-length datasets that you have kindly referenced: both $\\texttt{KSig-PDE}$ and $\\texttt{KSig-RFSF}$ failed to compile due to excessive memory requirements, and even the latter required significantly more memory than $\\texttt{PowerSig}$ and more than our current resources could handle. Our reproducing of prior results for comparison was further complicated by undocumented yet essential and nontrivial preprocessing steps and the absence of reference hyperparameter settings on the side of competing methods, making fair or even meaningful comparisons infeasible without extensive tuning and partial reimplementation. In particular, without such tuning or preprocessing, the generated signature kernel matrices, whether from near-exact approximations or low-rank variants alike, exhibited extremely high condition numbers on the given data, which rendered standard regression and classification tasks numerically unstable or entirely infeasible.\n\nThat said, we were already able to benchmark on the kernel regression dataset on bitcoin pricing showcased in [Salvi et al: \"The Signature Kernel is the Solution of a Goursat PDE\", SIAM J. Math. Data Sci., Vol. 3 (3)], where $\\texttt{PowerSig}$ has achieved increased accuracy (1.34% MAPE) compared to $\\texttt{KSigPDE}$ (1.41%) at a fraction of the memory, and we are actively working very hard to close this gap more convincingly as soon as possible and set up meaningful comparisons on selected UEA datasets and to also identify alternatives with tractable length and dimensionality where all competing methods can be run reliably. We aim to report updates on these benchmarks here as soon as possible, and certainly want to include at least a few such comparisons in the camera-ready version.\n\n   4. To make a compelling case for $\\texttt{PowerSig}$, I suggest that the authors provide numerical studies on at least a few learning tasks; e.g., Multivariate UEA Datasets, https://github.com/tgcsaba/KSig and compare against alternative approaches, Given the efficiency of $\\texttt{PowerSig}$, such studies should be feasible within a reasonable computational budget. This will also provide some insights into whether the benefits of (almost) exact signature kernels translate into superior predictive performance in practice.\n\nReply: We much appreciate this suggestion and have incorporated our response into the answer to the previous point.\n\n    5. Numerical studies on higher-dimensional paths (say d=10, 50) to confirm that the runtime grows linearly in d and accuracy remains stable.\n\nReply: Thank you for the suggestion. We have conducted the desired experiment on runtime vs time-series dimension on fixed-length (= 512) Brownian motion paths of increasing dimension, ranging from 2 to 4096, where we observed stable accuracy and confirmed the expected linear growth on runtime; indeed, we observed near perfect linear runtime growth from d=64 to d=8192, and sublinear growth outside.\n\n    6. Numerical studies showing how P should scale as a function of roughness and sequence length to maintain a target accuracy.\n\nReply: Thank you for suggesting this. Based on the discussions in Appendix B.3, we included in our original submission a utility function to address exactly this point: estimate_required_order() (in powersig.jax.algorithm.py) takes as input time series $\\(x_\\mu \\mid \\mu = 1,\\ldots, m)$, $(y_\\nu \\mid \\nu = 1,\\ldots, n)$, and a target error tolerance, and returns a recommended global truncation depth (capped at 64) for computing the cross-Gram matrix $(K_{x_\\mu,y_\\nu})\\in\\mathbb{R}^{m\\times n}$. Empirically, and in line with our theoretical discussions, the required truncation depth increases with both time-series roughness and sequence length. In practice, we observed a stable operational range between orders 7 and 15, beyond which further increases brought diminishing returns; instead, rescaling the time series was generally a more effective means to restore numerical stability. We will add a corresponding remark in the main text and explicitly refer to this utility function."}, {"review_id": "cKJeUkdvGX", "reviewer": "Reviewer_p9wB", "summary": "This paper presents PowerSig, a space-efficient approach to compute signature kernels. The signature kernel is a tool to analyze time series, which was previously generally computed via dynamic programming or as a solution to a Goursat PDE problem. Instead of solving a global PDE problem, this paper proposes to partition the unit square and solve many local PDE problems instead. These are reformulated as integral equations which admit local power series expansion solutions. Coefficients are first computed for the tile $T_{1,1}$ in the corner, after which boundary conditions are expanded alongside diagonals $T_{i, j: i+j=C}$ until the opposite corner tile is reached. The proposed method has similar time complexity compared to existing approaches, but requires significantly less memory especially for long sequences.", "questions": "1. Could the algorithm be modified to start at any of the four corner tiles? If not, why? If yes, could you simultaneously start from every corner and process tiles in parallel towards the middle of the square?\n\n2. Can you combine aspects of both your method and polysigkernel, or are they incompatible with each other?\n\n*Minor comments/typos*\n\nL56: \"solition\"\n\nSection 2.1: Why use a notation with 2 different \"l\" symbols ($\\ell, l$)?", "limitations": "Yes", "rating": 5, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 2, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**Strengths**\n\n- The proposed approach, due to its lower space complexity, allows to efficiently scale signature kernel computations to very long or rough time series, which is validated by numerical experiments.\n\n- The rapid convergence of the power series leads to highly accurate results, which may surpass existing PDE-based methods.\n\n- I find the overall approach elegant as it converts a complex global problem into simpler local ones.\n\n**Weaknesses**\n\n- For someone not already closely familiar with the signature kernel (in particular as a solution to a Goursat PDE), the mathematical derivations can be difficult to follow.\n\n- To further demonstrate the usefulness of the proposed method and reinforce the claim that it is suitable for large-scale applications in different domains, it could have been interesting to show how it could allow solving a previously intractable real-world problem.", "was_revised": true, "final_justification": "This submission introduces a space-efficient approach to compute signature kernels. In their rebuttal, the authors addressed the clarity issue I raised and started applying their method to more complex scenarios.", "rebuttal": "Thank you very much for your time and your very valuable feedback and suggestions. In what follows, we address each point separately, hopefully answering and clarifying all items. If anything further can be done from our side to facilitate this process, kindly let us know.\n\n1. For someone not already closely familiar with the signature kernel (in particular as a solution to a Goursat PDE), the mathematical derivations can be difficult to follow.\n\nReply: Thank you for pointing this out. We agree and have revised Section 2.1 accordingly by providing additional context and references on the Goursat PDE formulation of the signature kernel. We will also streamline the exposition of Sections 2.2 and 2.3 to improve readability, in particular by adding short \"road‑map\" paragraphs that clearly outline the purpose of each step before presenting the detailed technical derivations.\n\n2. To further demonstrate the usefulness of the proposed method and reinforce the claim that it is suitable for large-scale applications in different domains, it could have been interesting to show how it could allow solving a previously intractable real-world problem.\n\nReply: We agree. In response to this, we have benchmarked $\\texttt{PowerSig}$ on the kernel regression task for bitcoin pricing presented in [Salvi et al: \"The Signature Kernel is the Solution of a Goursat PDE\", SIAM J. Math. Data Sci., Vol. 3 (3)], where our $\\texttt{PowerSig}$ achieved improved accuracy over $\\texttt{KSig-PDE}$ while using only a fraction of the memory. This experiment provides an initial indication that $\\texttt{PowerSig}$ enables accuracy gains in practical kernel computations on real-world datasets that can become prohibitively large for existing signature kernel methods. We are currently working to demonstrating broader applicability in the spirit of your suggestion. So far, this effort has been hindered by severe resource limitations: competing methods such as $\\texttt{KSig-PDE}$ and $\\texttt{KSig-RFSF}$ failed to compile or run due to excessive memory requirements, and their default configurations proved unstable (e.g., generating ill-conditioned Gram matrices) without extensive, undocumented preprocessing and hyperparameter tuning.\n\nThat said, we are actively working to expand our benchmarks before the end of the rebuttal period to include additional real-world problems that have previously been out of reach for kernel-based methods. Our goal is to include at least a few such large-scale applications in the camera-ready version.\n\n3. Could the algorithm be modified to start at any of the four corner tiles? If not, why? If yes, could you simultaneously start from every corner and process tiles in parallel towards the middle of the square?\n\nReply: Thank you for the question and its follow-up. Indeed, the algorithm can only start at tiles where the boundary conditions of the Goursat PDE are fully specified, which in our formulation is exclusively the bottom-left corner tile. The reason is that propagating the signature kernel via the Goursat PDE on a given tile (say, $[\\sigma_k, \\sigma_{k+1}]\\times[\\tau_l,\\tau_{l+1}]$) requires knowledge of its initial conditions on both associated time-front edges ($\\{\\sigma_k\\}\\times[\\tau_l,\\tau_{l+1}]$ and $[\\sigma_k, \\sigma_{k+1}]\\times\\{\\tau_l\\}$), which intersect at the tile's bottom-left corner. These initial values are only known a priori on the global boundary edges $\\{0\\}\\times[0,1]$ and $[0,1]\\times\\{0\\}$, where the kernel is identically $1$ by definition. These edges intersect at the bottom-left tile of the grid, making it the only valid starting point under the current scheme. To initiate propagation from a different corner, one would need to perform a tile rotation by applying a suitable change-of-variables to the time series data. We will add a remark in our revised exposition to clarify this. \n\n4. Can you combine aspects of both your method and $\\texttt{polysigkernel}$, or are they incompatible with each other?\n\nReply: Thank you for this question. Both our method ($\\texttt{PowerSig}$) and the concurrent $\\texttt{polysigkernel}$ essentially target the same analytic object—the solution of the signature kernel's Goursat PDE—and yield discretizations via recursively propagated local power series expansions. They differ mainly conceptually and in their mathematical derivations: $\\texttt{polysigkernel}$ constructs its localised series expansions from an asymptotic representation involving integrated Bessel functions and analytic boundary conditions, while our approach develops and employs a direct functional-analytic perspective using a Neumann series representation based on the spectral properties of the integral operator that propagates the signature kernel via its characterising Goursat PDE; we elaborate on this distinction in a newly added comment in the manuscript.\n\nDespite this methodological distinction, the resulting algorithms are structurally very similar (though not entirely, and to our advantage, as newly added empirical benchmarks indicate), sharing recursion schemes, GPU-parallelizable wavefront propagation, and similarly rapid factorial convergence of the propagated power series. The implementations of our methods are indeed compatible: one could, for instance, replace the Neumann update in our solver with $\\texttt{polysigkernel}$'s Chebyshev step on selected tiles, or incorporate our low‑memory propagation update into their code. That said, given the close alignment, we believe that such combinations would offer little additional benefit. Motivated by your question, however, we see greater potential in selectively integrating finite-difference or FEM-based local solvers of the Goursat PDE (as used in $\\texttt{polysigkernel}$) on particularly challenging subclusters of the tile grid. This hybrid approach could enhance robustness while maintaining memory efficiency and parallel scalability. We will include a brief discussion of this in the revised manuscript."}, {"review_id": "WULv3WOvSx", "reviewer": "Reviewer_4Qtm", "summary": "This paper introduces PowerSig, a new method for computing the signature kernel—a powerful but computationally expensive kernel function for comparing time series based on rough path theory. The key idea is to avoid the traditional heavy global computation (e.g., solving a dense PDE or using dynamic programming) and instead break the kernel domain into small tiles, then compute local power series approximations within each tile using Neumann series expansions. Each tile’s kernel value is computed recursively using only its immediate neighbors, which significantly reduces memory usage and makes the method highly parallelizable. Because the power series converge quickly and can be adaptively truncated, the method maintains high accuracy even on very long or rough time series. The authors demonstrate strong empirical improvements in speed, memory efficiency, and numerical stability over existing state-of-the-art solvers.", "questions": "Do you have any theoretical insights or empirical heuristics that guide truncation depth per tile (e.g., based on local values of \\rho)?", "limitations": "While the results for low-Hurst fractional Brownian motion are impressive, I wonder about more pathological edge cases. Are there regimes (e.g., extremely high-dimensional paths or adversarial oscillations in x and y) where the power series becomes unstable or requires very high truncation orders?", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "I found the technical contribution practical. The use of localized Neumann series to approximate the signature kernel is a smart way to circumvent the global memory bottleneck, and the paper demonstrates that this leads to substantial improvements in scalability. The method is grounded in solid math (with rigorous derivations and proofs), and the experimental results back up the claims well -- especially the ability to handle time series with over a million points on a single GPU. I also appreciate how the algorithm leverages the geometry of the time series in a structured, recursive way. This feels like a well-engineered solution to a real bottleneck in signature kernel applications. One area I wish the paper had gone a bit deeper is in the comparison and discussion with related or emerging alternatives. For example, the concurrent work by Cass et al. is mentioned but not actually benchmarked or contrasted. I feel like this is a critical piece for the completeness of this work. It would help contextualize PowerSig’s practical impact more clearly if the authors could situate it relative to other fast approximations or sparse kernel techniques (beyond KSig). Also, while the math is clean, some of the technical exposition -- particularly in Section 2.2 -- feels heavy and could benefit from more intuitive scaffolding or diagrams to make the recursion scheme easier to follow on a first read.", "was_revised": true, "final_justification": "The authors have successfully addresses my major concerns, i.e. benchmarking with a similar methods, providing more guidance on how to choose truncation depth, etc. I will be happy to see all these being incorporated with the camera-ready version, which makes the paper more compelling and comprehensive.", "rebuttal": "Thank you very much for your time and the helpful comments and questions. In what follows, we address each point separately, hopefully answering any concerns you may have. Kindly let us know if you wish any further clarifications or answers.\n\n1. One area I wish the paper had gone a bit deeper is in the comparison and discussion with related or emerging alternatives. For example, the concurrent work by Cass et al. is mentioned but not actually benchmarked or contrasted. I feel like this is a critical piece for the completeness of this work.\n\nReply: We have since benchmarked and contrasted our method against the concurrent work of Cass et al. Both approaches were developed independently and released within two weeks of each other, with both preprints carrying the same arXiv month stamp. To address the understandable request for a direct comparison of our methods, our revised manuscript now contains a new subsection entitled \"Benchmark against polysigkernel\", accompanied by an additional figure panel showcasing performance on joint benchmarks. Using the publicly available JAX implementation of $\\texttt{polysigkernel}$, we ran identical experiments on computation time and memory usage over two standard benchmarks: $(\\alpha)$ 2D Brownian motion paths of length 2 to 512, and $(\\beta)$ fixed-length (= 512) paths of increasing dimension, ranging from 2 to 4096. These experiments, which we are happy to detail further upon request and will include with the camera-ready version, show that both solvers achieve essentially identical accuracy throughout while $\\texttt{PowerSig}$ performs better in runtime and memory usage. In benchmark $(\\alpha)$, $\\texttt{PowerSig}$ was on average 3.1% faster and required 23% less memory than $\\texttt{polysigkernel}$ on oscillatory time series of length 512 (averaged over 10 iid samples). In benchmark $(\\beta)$, $\\texttt{PowerSig}$ exhibited a comparable computational advantage, being 8.5% faster on average at dimension 512, with similar memory usage.\n\nOur added subsection also clarifies why these empirical differences occur and how our underlying methods differ in concept despite propagating the same local power‑series map: $\\texttt{Polysigkernel}$ starts from an explicit Riemann‑function formula to rewrite the Goursat solution as a sum of integrated modified Bessel functions $I_{0}$ [Cass et al., Thm. 3.1 (Polyanin)] which, after expanding each $I_{0}$ in its well‑known power series, allows the kernel on every rectangle to be reduced (using integration‑by‑parts and a change of variables) to two univariate centred power series in $s$ and $t$, for which the resulting coefficient arrays can then be generated by the closed‑form recurrence $\\Lambda_{C;[\\sigma_{k},\\sigma_{k+1}]\\times[\\tau_{l},\\tau_{l+1}]}$ (their Props. 3.1, A.1, A.2). Truncating these expansions up to some fixed degree $N$ yields a recursive polynomial approximation scheme ($\\texttt{polysigkernel}$) for which the authors provide local and global truncation error bounds (their Prop. 3.4 and Thm. 3.2). A lower-performant Chebyshev‑interpolation scheme for polynomial boundary data is also discussed in their Sect. 3.3. Our algorithm, by contrast, approaches the same PDE from a functional-analytic view: Using the Volterra formulation of the Goursat problem, we express the kernel propagator through an integral operator with favourable spectrum (Lemma 2.3) that permits a concise three-step Neumann recursion (Prop. 2.5) whose iterates build the same local power series coefficient-by-coefficient in a way that allows for their immediate extraction on the fly (Prop. 2.8). This Neumann expansion terminates as soon as the next term is below machine precision, so iteration depth that can vary per tile and no global cut‑off $N$ is required; in practice 5–8 iterations suffice except on the roughest tiles, resulting in a memory footprint that stores only one tile strip and is strictly lower than the quadratic complexity held by the given $\\texttt{polysigkernel}$ implementation.\n\nIn summary, $\\texttt{polysigkernel}$ obtains each tile's coefficients from a closed‑form Bessel (or Chebyshev) recurrence evaluated up to a user‑chosen degree $N$, whereas $\\texttt{PowerSig}$ builds the same coefficients adaptively via a Neumann iteration that stops once the next term drops below machine precision. Because both procedures implement exactly the same propagation map, they are analytically interchangeable on every tile.  Their distinct constructions, however, naturally lead to different implementations and practical profiles, and our head‑to‑head benchmarks show a modest but consistent performance edge for $\\texttt{PowerSig}$.\n\n2. It would help contextualize $\\texttt{PowerSig}$'s practical impact more clearly if the authors could situate it relative to other fast approximations or sparse kernel techniques (beyond KSig).\n\nReply: Thank you for this suggestion. We agree and are currently working to implement it by benchmarking $\\texttt{PowerSig}$ against the recent specialised low-rank and RFF-based approaches introduced in [Tóth, Oberhauser, and Szabó: \"Random Fourier Signature Features\", SIAM J. Math. Data Sci., Vol. 7 (1)]. We expect to have results available in the coming days and will include them in the camera-ready version. Any additional suggestions of specific sparse kernel techniques or alternative fast approximations suitable for comparison are welcome.\n\n3. Also, while the math is clean, some of the technical exposition -- particularly in Section 2.2 -- feels heavy and could benefit from more intuitive scaffolding or diagrams to make the recursion scheme easier to follow on a first read.\n\nReply: Thank you for pointing this out. We share your concern and will address this in two ways in the camera-ready version: a) the summarising schematic Figure 3 (currently in Appendix B), which illustrates our recursion scheme and its information flow, will be moved to the main text and placed alongside the current summary paragraph at the start of Section 2.2. We will also expand the figure caption to walk the reader through a full propagation step, clarifying the role of each object before any algebra appears. b) Additionally, we will add a short introductory paragraph (four to five sentences) that informally explains why only the south‑ and west‑edge data of a tile enter the update, thereby linking the figure more explicitly to the characteristic structure of the Goursat PDE. With these adjustments, the reader will first encounter an intuitive scaffold before engaging with the formal derivation, while all mathematical information remains unchanged.\n\n4. Do you have any theoretical insights or empirical heuristics that guide truncation depth per tile (e.g., based on local values of $\\rho$)?\n\nReply: Yes, we do. Proposition 2.8 (cf. also Appendix B.3) shows that the approximation error of our expansion on a tile of side-length $\\Delta$ at truncation depth $N$ is bounded by $\\sum_{i,j \\,:\\, \\max(i,j)\\geq N+1}a_{ij}$ for $a_{ij}:=(|\\rho|\\Delta^2)^{\\min(i,j)}/(\\min(i,j)!)^2$, where $\\rho$ is the local value of the Goursat PDE coefficient function on that tile. As described in the appendix, we have taken this inequality as the basis for a utility function, implemented as estimate_required_order() in the provided file powersig.jax.algorithm.py, which takes as input time series $(x_\\mu \\mid \\mu = 1,\\ldots, m)$ and $(y_\\nu \\mid \\nu = 1,\\ldots, n)$, along with a desired error tolerance, and suggests an optimal global truncation depth (capped at 64) for computing the cross-Gram matrix $(K_{x_\\mu,y_\\nu})\\in\\mathbb{R}^{m\\times n}$. In our use cases, we found the numerical sweet spot to lie between truncation orders 7 and 15; beyond this range, rescaling the time series, rather than increased truncation depth, was the more effective way to restore stability.\n\n5. While the results for low-Hurst fractional Brownian motion are impressive, I wonder about more pathological edge cases. Are there regimes (e.g., extremely high-dimensional paths or adversarial oscillations in x and y) where the power series becomes unstable or requires very high truncation orders? \n\nReply: This points to an important aspect, thank you for raising it. As discussed in Appendix B.3, the numerical efficiency and stability of our tile-based approximation are closely tied to the behaviour of the inner product coefficient $\\rho_{x,y}(s,t) =\\langle \\frac{dx}{ds},\\frac{dy}{dt} \\rangle$ that governs the local Goursat PDE. Instabilities can occur when $\\rho_{x,y}$ exhibits large spikes on isolated tiles, pushing local truncation orders high. This typically arises when the time series have strongly heterogeneous increments, e.g. a few exceptionally large differences amid otherwise small ones, leading to extreme variations in $|\\rho_{x,y}|$ across tiles. Such behaviour can break the 'generic situation' that local power series decay rapidly and uniformly.\n\nAdversarial oscillations (in the sense of irregular or bursty increments) can trigger this, particularly when isolated increments dominate in magnitude. High-dimensional paths, by contrast, are usually benign unless they exhibit the same large-increment structure, since the dimensionality enters only through $\\rho_{x,y}$. In practice, these edge cases can be flagged when our internal utility function falls outside the typical truncation range of $[7,15]$. We have found that stability is then often recoverable via simple, invertible preprocessing (e.g., normalization or whitening) to compress the range of $|\\rho_{x,y}|$ into a numerically stable regime."}]

Author rebuttal:

Meta-review: This manuscript proposes a method for efficiently working with the signature kernel, particularly when the underlying time-series data is long or high-dimensional. A key technical innovation is the use of local power series expansions (with appropriate surrounding machinery). 

The reviewers all highlight the efficacy of the scheme; it is interesting and clearly effective. The numerical experiments clearly highlight this. Moreover, in the manuscript it is accompanied by theory well matched to its practical use.

Some minor weaknesses are in the presentation (e.g., being quite dependent on a lot of "prior" information) and, initially, the scope of some experiments. The latter was addressed during the rebuttal; the former may be somewhat intrinsic.
Decision: Accept (poster)